library("knitr") # for knitting stuff
library("lme4") # for linear mixed effects models
library("broom") # for tidying model results
library("brms") # Bayesian regression with Stan
library("corrr") # for tidy correlation matrix
library("xtable") # for latex tables
library("jsonlite") # for reading json files
library("janitor") # cleans stuff
library("ggtern") # for ternary plots
library("Hmisc") # bootstrapped confidence intervals
library("tidybayes") # for Bayesian data analysis
library("png") # adding pngs to images
library("grid") # functions for dealing with images
library("egg") # for geom_custom
library("tidyverse") # for wrangling, plotting, etc. # root mean squared error
rmse = function(x, y){
return(sqrt(mean((x - y)^2)))
}
actual_counterfactual_plot = function(x, ylabel) {
p = ggplot(data = x,
mapping = aes(x = clipindex,
y = value,
fill = colorindex,
group = model)) +
stat_summary(geom = "bar",
fun.y = "mean",
color = "black",
position = position_dodge(0.7),
width = 0.7) +
stat_summary(geom = "linerange",
fun.data = "mean_cl_boot",
position = position_dodge(0.7),
size = 1) +
geom_point(data = x %>% mutate(value = ifelse(model == "rating", value, NA)),
position = position_jitterdodge(dodge.width = 0.7, jitter.width = 0.2),
shape = 21,
color = "black",
size = 1,
alpha = 0.2) +
scale_fill_manual(values = c("red2", "green2", "black")) +
facet_grid(index_actual ~ index_counterfactual) +
geom_text(data = df.text %>% filter(!is.na(label)),
mapping = aes(x = x, y = y, label = as.character(label)),
size = 6,
position = position_dodge(width = .7),
hjust = 0.5) +
coord_cartesian(xlim = c(0.5, 2.5),
ylim = c(-5, 105),
expand = F,
clip = "off") +
labs(y = ylabel) +
theme_bw() +
theme(
text = element_text(size = 20),
legend.position = "none",
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(margin = margin(0, 10, 0, 0)),
panel.spacing.y = unit(c(.8), "cm"),
plot.margin = unit(c(0.2, 0.2, .8, .2), "cm"),
axis.title.x = element_blank(),
panel.grid = element_blank(),
strip.background = element_blank(),
panel.border = element_rect(colour = "black")
)
print(p)
}
actual_counterfactual_threeballs_plot = function(x, ylabel) {
p = ggplot(data = x,
mapping = aes(x = ball,
y = value,
group = model,
fill = colorindex)) +
stat_summary(fun.y = "mean",
geom = "bar",
color = "black",
position = position_dodge(0.9),
width = 0.9) +
stat_summary(fun.data = "mean_cl_boot",
geom = "linerange",
color = "black",
position = position_dodge(0.9)) +
geom_point(data = x %>% mutate(value = ifelse(model == "rating", value, NA)),
position = position_jitterdodge(jitter.width = 0.3, jitter.height = 0, dodge.width = 0.9),
shape = 21,
color = "black",
size = 1,
alpha = 0.15) +
facet_wrap(~clip, ncol = 8) +
geom_text(data = df.text,
mapping = aes(x = x,
y = y,
label = as.character(label)),
size = 5,
position = position_dodge(width = .9),
hjust = 0.5,
fontface = "bold") +
coord_cartesian(xlim = c(0.5, 2.5),
ylim = c(-5, 120),
expand = F,
clip = "off") +
scale_y_continuous(breaks = seq(0, 100, 25)) +
scale_fill_manual(values = c("red2", "green2", "black")) +
scale_color_manual(values = c("red2", "green2", "black")) +
labs(y = ylabel) +
theme_bw() +
theme(
text = element_text(size = 16),
legend.position = "none",
axis.title.x = element_blank(),
panel.spacing.y = unit(c(.3), "cm"),
panel.grid = element_blank(),
strip.background = element_blank(),
strip.text = element_blank(),
panel.border = element_rect(colour = "black")
)
print(p)
}# causal judgments
df.exp1.causal = read.csv("../../data/experiment1_causal.csv",
header = T,
stringsAsFactors = F)
# counterfactual judgments
df.exp1.counterfactual = read.csv("../../data/experiment1_counterfactual.csv",
header = T,
stringsAsFactors = F)
# Information about each clip
df.exp1.clipinfo = tibble(
clip = 1:18,
outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1),
outcome_counterfactual = c(0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1),
index_actual = rep(c("actual miss",
"actual close",
"actual hit"), each = 6),
index_counterfactual = rep(rep(c("counterfactual miss",
"counterfactual close",
"counterfactual hit"), each = 2), 3)
) %>%
mutate(index_actual = factor(index_actual, levels = c("actual miss",
"actual close",
"actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")))
# Structure data
df.exp1.causal.long = df.exp1.causal$data %>%
strsplit("\n") %>%
unlist() %>%
strsplit(",") %>%
unlist() %>%
matrix(ncol = 36, byrow = T) %>%
as.data.frame(stringsAsFactors = F) %>%
mutate_all("as.numeric") %>%
setNames(str_c(c("clip_", "rating_"), rep(1:18, each = 2))) %>%
mutate(participant = 1:nrow(.)) %>%
gather("index", "value", -participant) %>%
separate(index, into = c("index", "clip")) %>%
spread(index, value) %>%
left_join(df.exp1.clipinfo,
by = "clip") %>%
arrange(participant, clip) %>%
ungroup()
df.exp1.counterfactual.long = df.exp1.counterfactual$data %>%
strsplit("\n") %>%
unlist() %>%
strsplit(",") %>%
unlist() %>%
matrix(ncol = 36, byrow = T) %>%
as.data.frame(stringsAsFactors = F) %>%
mutate_all("as.numeric") %>%
setNames(str_c(c("clip_", "rating_"), rep(1:18, each = 2))) %>%
mutate(participant = 1:nrow(.)) %>%
gather("index", "value", -participant) %>%
separate(index, into = c("index", "clip")) %>%
spread(index, value) %>%
left_join(df.exp1.clipinfo,
by = "clip") %>%
arrange(participant, clip) %>%
ungroup()# counterfactual condition
df.exp1.counterfactual %>%
separate(extra,
into = c("gender", "age", "difficulty", "smoothness", "time"),
sep = ",") %>%
mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>%
select(gender:time) %>%
mutate_if(is.character, ~ as.numeric(.)) %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2))
# causal condition
df.exp1.causal %>%
separate(extra,
into = c("gender", "age", "difficulty", "smoothness", "time"),
sep = ",") %>%
mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>%
select(gender:time) %>%
mutate_if(is.character, ~ as.numeric(.)) %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) n age_mean age_sd n_female time_mean time_sd
1 41 32 8.8 19 7.5 3.12
n age_mean age_sd n_female time_mean time_sd
1 41 33 10.3 21 9.65 6.58
# read model predictions
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "2ball")]
df.exp1.model = files %>%
map(~ read.csv(str_c(path, .))) %>%
reduce(rbind) %>%
set_names(c("clip", "noise", "prediction")) %>%
left_join(df.exp1.clipinfo, by = "clip") %>%
mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))
# calculate mean counterfactual judgments
df.exp1.counterfactual.means = df.exp1.counterfactual.long %>%
group_by(clip) %>%
summarize(rating_mean = mean(rating),
rating_low = smean.cl.boot(rating)[2],
rating_high = smean.cl.boot(rating)[3]) %>%
ungroup() %>%
left_join(df.exp1.clipinfo, by = "clip")
# find noisy simulation model that correlates best with the mean counterfactual judgments
df.exp1.counterfactual.model = df.exp1.model %>%
group_by(noise) %>%
nest() %>%
mutate(correlation = map(data, ~ cor(.x$prediction,
df.exp1.counterfactual.means$rating_mean)),
correlation = unlist(correlation)) %>%
ungroup() %>%
filter(correlation == max(correlation)) %>%
unnest(data) %>%
mutate(prediction = prediction * 100)df.exp1.causal.model = df.exp1.counterfactual.long %>%
group_by(clip) %>%
summarize(empirical = mean(rating)) %>%
ungroup() %>%
left_join(df.exp1.counterfactual.model %>%
select(-c(noise, correlation)),
by = "clip") %>%
rename(simulation = prediction) %>%
mutate(empirical = ifelse(outcome_actual == 1, 100 - empirical, empirical),
simulation = ifelse(outcome_actual == 1, 100 - simulation, simulation))df.plot = df.exp1.counterfactual.long %>%
mutate(rating = abs(rating),
clipindex = rep(c(1, 2), nrow(.)/2)) %>%
left_join(df.exp1.counterfactual.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
gather("model", "value", c(rating, prediction)) %>%
mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
colorindex = ifelse(model != "rating", 3, colorindex),
colorindex = as.factor(colorindex),
index_actual = factor(index_actual,levels = c("actual miss", "actual close", "actual hit")),
index_counterfactual = factor(index_counterfactual,levels = c("counterfactual miss", "counterfactual close", "counterfactual hit")),
model = factor(model, levels = c("rating", "prediction"))) %>%
mutate(clipindex = ifelse(clip == 11, 2, clipindex),
clipindex = ifelse(clip == 12, 1, clipindex)) #swap clips 11 and 12
df.text = df.plot %>%
expand(index_actual, index_counterfactual, clipindex, model) %>%
mutate(
label = rep(1:18, each = 2),
label = ifelse(model != "rating", NA, label),
y = -15,
x = clipindex,
colorindex = NA
)
actual_counterfactual_plot(df.plot, "counterfactual judgment")Warning: Removed 738 rows containing missing values (geom_point).
df.exp1.counterfactual.means %>%
left_join(df.exp1.counterfactual.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
summarize(simulation_r = cor(rating_mean, prediction),
simulation_rmse = rmse(rating_mean, prediction)) %>%
mutate_all(~ round(., 2))# A tibble: 1 x 2
simulation_r simulation_rmse
<dbl> <dbl>
1 0.98 12.0
model.name = c("empirical")
# model.name = c("simulation")
df.plot = df.exp1.causal.long %>%
mutate(
rating = abs(rating),
clipindex = rep(c(1, 2), nrow(.) / 2)
) %>%
left_join(df.exp1.causal.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
gather("model", "value", c(rating, simulation, empirical)) %>%
filter(model %in% c(model.name, "rating")) %>%
mutate(model = factor(model, levels = c("rating", model.name))) %>%
mutate(
colorindex = ifelse(outcome_actual == 1, 2, 1),
colorindex = ifelse(model != "rating", 3, colorindex),
colorindex = as.factor(colorindex),
index.actual = factor(index_actual,
levels = c("actual miss", "actual close", "actual hit")),
index.counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit"))
) %>%
mutate(
clipindex = ifelse(clip == 11, 2, clipindex),
clipindex = ifelse(clip == 12, 1, clipindex)
) # swap clips 11 and 12
df.text = df.plot %>%
expand(index_actual, index_counterfactual, clipindex, model) %>%
mutate(
label = rep(1:18, each = 2),
label = ifelse(model != "rating", NA, label),
y = -15,
x = clipindex,
colorindex = NA
)
actual_counterfactual_plot(df.plot, "causal judgment")Warning: Removed 738 rows containing missing values (geom_point).
df.exp1.causal.long %>%
mutate(rating = abs(rating)) %>%
group_by(clip,
outcome_actual,
outcome_counterfactual,
index_actual,
index_counterfactual) %>%
summarize(rating = mean(rating)) %>%
ungroup() %>%
left_join(df.exp1.causal.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
summarize(empirical_r = cor(rating, empirical),
empirical_rmse = rmse(rating, empirical),
simulation_r = cor(rating, simulation),
simulation_rmse = rmse(rating, simulation)) %>%
mutate_all(~ round(., 2))# A tibble: 1 x 4
empirical_r empirical_rmse simulation_r simulation_rmse
<dbl> <dbl> <dbl> <dbl>
1 0.96 8.57 0.94 16.8
# random intercepts model
fit.brm_exp1_causal = brm(formula = rating ~ 1 + index_actual + index_counterfactual + outcome_actual + (1 + index_actual + index_counterfactual + outcome_actual | participant),
data = df.exp1.causal.long,
cores = 2,
seed = 1,
file = "cache/brm_exp1_causal")
fit.brm_exp1_causal %>%
summary() Family: gaussian
Links: mu = identity; sigma = identity
Formula: rating ~ 1 + index_actual + index_counterfactual + outcome_actual + (1 + index_actual + index_counterfactual + outcome_actual | participant)
Data: df.exp1.causal.long (Number of observations: 738)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~participant (Number of levels: 41)
Estimate
sd(Intercept) 13.95
sd(index_actualactualclose) 3.40
sd(index_actualactualhit) 5.51
sd(index_counterfactualcounterfactualclose) 11.02
sd(index_counterfactualcounterfactualhit) 21.28
sd(outcome_actual) 22.11
cor(Intercept,index_actualactualclose) -0.18
cor(Intercept,index_actualactualhit) -0.03
cor(index_actualactualclose,index_actualactualhit) 0.11
cor(Intercept,index_counterfactualcounterfactualclose) 0.15
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 0.03
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 0.03
cor(Intercept,index_counterfactualcounterfactualhit) -0.43
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 0.10
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 0.17
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) -0.05
cor(Intercept,outcome_actual) -0.48
cor(index_actualactualclose,outcome_actual) -0.02
cor(index_actualactualhit,outcome_actual) -0.26
cor(index_counterfactualcounterfactualclose,outcome_actual) -0.17
cor(index_counterfactualcounterfactualhit,outcome_actual) -0.44
Est.Error
sd(Intercept) 3.12
sd(index_actualactualclose) 2.59
sd(index_actualactualhit) 4.10
sd(index_counterfactualcounterfactualclose) 4.46
sd(index_counterfactualcounterfactualhit) 3.91
sd(outcome_actual) 3.93
cor(Intercept,index_actualactualclose) 0.37
cor(Intercept,index_actualactualhit) 0.36
cor(index_actualactualclose,index_actualactualhit) 0.38
cor(Intercept,index_counterfactualcounterfactualclose) 0.29
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 0.37
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 0.36
cor(Intercept,index_counterfactualcounterfactualhit) 0.19
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 0.36
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 0.34
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 0.28
cor(Intercept,outcome_actual) 0.18
cor(index_actualactualclose,outcome_actual) 0.36
cor(index_actualactualhit,outcome_actual) 0.36
cor(index_counterfactualcounterfactualclose,outcome_actual) 0.28
cor(index_counterfactualcounterfactualhit,outcome_actual) 0.19
l-95% CI
sd(Intercept) 8.14
sd(index_actualactualclose) 0.14
sd(index_actualactualhit) 0.18
sd(index_counterfactualcounterfactualclose) 2.13
sd(index_counterfactualcounterfactualhit) 14.23
sd(outcome_actual) 15.17
cor(Intercept,index_actualactualclose) -0.79
cor(Intercept,index_actualactualhit) -0.70
cor(index_actualactualclose,index_actualactualhit) -0.64
cor(Intercept,index_counterfactualcounterfactualclose) -0.41
cor(index_actualactualclose,index_counterfactualcounterfactualclose) -0.70
cor(index_actualactualhit,index_counterfactualcounterfactualclose) -0.67
cor(Intercept,index_counterfactualcounterfactualhit) -0.75
cor(index_actualactualclose,index_counterfactualcounterfactualhit) -0.63
cor(index_actualactualhit,index_counterfactualcounterfactualhit) -0.52
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) -0.60
cor(Intercept,outcome_actual) -0.78
cor(index_actualactualclose,outcome_actual) -0.68
cor(index_actualactualhit,outcome_actual) -0.84
cor(index_counterfactualcounterfactualclose,outcome_actual) -0.66
cor(index_counterfactualcounterfactualhit,outcome_actual) -0.76
u-95% CI
sd(Intercept) 20.41
sd(index_actualactualclose) 9.69
sd(index_actualactualhit) 15.38
sd(index_counterfactualcounterfactualclose) 19.72
sd(index_counterfactualcounterfactualhit) 29.22
sd(outcome_actual) 30.76
cor(Intercept,index_actualactualclose) 0.59
cor(Intercept,index_actualactualhit) 0.65
cor(index_actualactualclose,index_actualactualhit) 0.78
cor(Intercept,index_counterfactualcounterfactualclose) 0.70
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 0.70
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 0.70
cor(Intercept,index_counterfactualcounterfactualhit) -0.02
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 0.73
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 0.77
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 0.47
cor(Intercept,outcome_actual) -0.08
cor(index_actualactualclose,outcome_actual) 0.68
cor(index_actualactualhit,outcome_actual) 0.50
cor(index_counterfactualcounterfactualclose,outcome_actual) 0.41
cor(index_counterfactualcounterfactualhit,outcome_actual) -0.06
Rhat
sd(Intercept) 1.00
sd(index_actualactualclose) 1.00
sd(index_actualactualhit) 1.00
sd(index_counterfactualcounterfactualclose) 1.00
sd(index_counterfactualcounterfactualhit) 1.00
sd(outcome_actual) 1.01
cor(Intercept,index_actualactualclose) 1.00
cor(Intercept,index_actualactualhit) 1.00
cor(index_actualactualclose,index_actualactualhit) 1.00
cor(Intercept,index_counterfactualcounterfactualclose) 1.00
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 1.00
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 1.00
cor(Intercept,index_counterfactualcounterfactualhit) 1.00
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 1.01
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 1.01
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 1.00
cor(Intercept,outcome_actual) 1.01
cor(index_actualactualclose,outcome_actual) 1.02
cor(index_actualactualhit,outcome_actual) 1.01
cor(index_counterfactualcounterfactualclose,outcome_actual) 1.00
cor(index_counterfactualcounterfactualhit,outcome_actual) 1.00
Bulk_ESS
sd(Intercept) 1599
sd(index_actualactualclose) 1553
sd(index_actualactualhit) 1127
sd(index_counterfactualcounterfactualclose) 805
sd(index_counterfactualcounterfactualhit) 1162
sd(outcome_actual) 1065
cor(Intercept,index_actualactualclose) 3783
cor(Intercept,index_actualactualhit) 4114
cor(index_actualactualclose,index_actualactualhit) 1953
cor(Intercept,index_counterfactualcounterfactualclose) 1567
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 828
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 1345
cor(Intercept,index_counterfactualcounterfactualhit) 1034
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 379
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 497
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 740
cor(Intercept,outcome_actual) 805
cor(index_actualactualclose,outcome_actual) 360
cor(index_actualactualhit,outcome_actual) 432
cor(index_counterfactualcounterfactualclose,outcome_actual) 900
cor(index_counterfactualcounterfactualhit,outcome_actual) 1742
Tail_ESS
sd(Intercept) 1756
sd(index_actualactualclose) 1582
sd(index_actualactualhit) 1691
sd(index_counterfactualcounterfactualclose) 1162
sd(index_counterfactualcounterfactualhit) 2131
sd(outcome_actual) 1768
cor(Intercept,index_actualactualclose) 3088
cor(Intercept,index_actualactualhit) 2881
cor(index_actualactualclose,index_actualactualhit) 2573
cor(Intercept,index_counterfactualcounterfactualclose) 2576
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 2032
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 2298
cor(Intercept,index_counterfactualcounterfactualhit) 2121
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 897
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 1619
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 1765
cor(Intercept,outcome_actual) 1632
cor(index_actualactualclose,outcome_actual) 810
cor(index_actualactualhit,outcome_actual) 1213
cor(index_counterfactualcounterfactualclose,outcome_actual) 1574
cor(index_counterfactualcounterfactualhit,outcome_actual) 2685
Population-Level Effects:
Estimate Est.Error l-95% CI
Intercept -18.33 3.36 -24.86
index_actualactualclose 4.31 3.40 -2.54
index_actualactualhit 4.06 4.84 -5.51
index_counterfactualcounterfactualclose -17.52 3.37 -24.36
index_counterfactualcounterfactualhit -61.37 4.38 -69.95
outcome_actual 92.36 5.11 82.21
u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -11.66 1.00 2026 2429
index_actualactualclose 10.94 1.00 2989 2747
index_actualactualhit 13.27 1.00 2531 2706
index_counterfactualcounterfactualclose -10.84 1.00 2898 2310
index_counterfactualcounterfactualhit -52.89 1.00 2210 2724
outcome_actual 102.38 1.00 2179 2885
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 30.09 0.89 28.42 31.83 1.00 3208 3142
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
df.exp1.causal.posterior = df.exp1.causal.long %>%
group_by(clip, index_actual, index_counterfactual, outcome_actual) %>%
summarize(value = mean(rating)) %>%
ungroup() %>%
add_fitted_draws(newdata = .,
model = fit.brm_exp1_causal,
re_formula = NA) %>%
ungroup() %>%
select(clip, .value, .draw) %>%
spread(clip, .value)
func_posterior_difference = function(data, clip1, clip2){
data %>%
mutate(difference = .data[[clip1]] - .data[[clip2]]) %>%
pull(difference) %>%
mean_hdci() %>%
mutate_if(is.numeric, ~round(., 2)) %>%
summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])")) %>%
print()
}
func_posterior_difference(data = df.exp1.causal.posterior,
clip1 = "8",
clip2 = "13") difference
1 (0.26 [-6.51, 7.31])
fit.brm_exp1_causal %>%
tidy() %>%
filter(str_detect(term, "b_")) %>%
mutate(term = str_remove(term, "b_"),
term = tolower(term)) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
rename(`lower 95% HDI` = lower,
`upper 95% HDI` = upper) %>%
xtable() %>%
print(include.rownames = F,
booktabs = T)% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:05:19 2019
\begin{table}[ht]
\centering
\begin{tabular}{lrrrr}
\toprule
term & estimate & std.error & lower 95\% HDI & upper 95\% HDI \\
\midrule
intercept & -18.33 & 3.36 & -23.81 & -12.78 \\
index\_actualactualclose & 4.31 & 3.40 & -1.40 & 9.83 \\
index\_actualactualhit & 4.06 & 4.84 & -4.01 & 11.88 \\
index\_counterfactualcounterfactualclose & -17.52 & 3.37 & -22.96 & -12.01 \\
index\_counterfactualcounterfactualhit & -61.37 & 4.38 & -68.54 & -54.30 \\
outcome\_actual & 92.36 & 5.11 & 84.04 & 100.96 \\
\bottomrule
\end{tabular}
\end{table}
df.exp2.causal_first = read.csv("../../data/experiment2_causal_first.csv", header = T, stringsAsFactors = F)
df.exp2.counterfactual_first = read.csv("../../data/experiment2_counterfactual_first.csv", header = T, stringsAsFactors = F)
# Information about each clip
df.exp2.clipinfo = tibble(
clip = 1:20,
outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1),
outcome_counterfactual = c(0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1),
index_actual = c(rep(c("actual miss", "actual close", "actual hit"), each = 6), rep("practice", 2)),
index_counterfactual = c(rep(rep(c("counterfactual miss", "counterfactual close", "counterfactual hit"), each = 2), 3), rep("practice", 2))
) %>%
mutate(index_actual = factor(index_actual, levels = c("actual miss",
"actual close",
"actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")))
# Structure data
df.exp2.causal_first.long = df.exp2.causal_first$data %>%
strsplit("\n") %>%
unlist() %>%
strsplit(",") %>%
unlist() %>%
matrix(ncol = 80, byrow = T) %>%
as.data.frame(stringsAsFactors = F) %>%
mutate_all("as.numeric") %>%
setNames(str_c(c("clip_", "rating_"), rep(1:40, each = 2))) %>%
mutate(participant = 1:nrow(.)) %>%
gather("index", "value", -participant) %>%
separate(index, into = c("index", "trial")) %>%
mutate(trial = as.numeric(trial)) %>%
spread(index, value) %>%
arrange(participant, trial) %>%
mutate(
question = rep(rep(c("causal", "counterfactual"), each = 20), max(participant)),
condition = "causal_first"
) %>%
left_join(df.exp2.clipinfo,
by = "clip") %>%
arrange(participant, question, clip) %>%
select(participant, question, clip, rating, everything())
df.exp2.counterfactual_first.long = df.exp2.counterfactual_first$data %>%
strsplit("\n") %>%
unlist() %>%
strsplit(",") %>%
unlist() %>%
matrix(ncol = 80, byrow = T) %>%
as.data.frame(stringsAsFactors = F) %>%
mutate_all("as.numeric") %>%
setNames(str_c(c("clip_", "rating_"), rep(1:40, each = 2))) %>%
mutate(participant = 1:nrow(.)) %>%
gather("index", "value", -participant) %>%
separate(index, into = c("index", "trial")) %>%
mutate(trial = as.numeric(trial)) %>%
spread(index, value) %>%
arrange(participant, trial) %>%
mutate(
question = rep(rep(c("counterfactual", "causal"), each = 20), max(participant)),
condition = "counterfactual_first"
) %>%
left_join(df.exp2.clipinfo,
by = "clip") %>%
arrange(participant, question, clip) %>%
select(participant, question, clip, rating, everything())
# combine data from both conditions
df.exp2.combined.long = df.exp2.causal_first.long %>%
bind_rows(df.exp2.counterfactual_first.long %>%
mutate(participant = participant + max(df.exp2.causal_first.long$participant)))# counterfactual condition
df.exp2.counterfactual_first %>%
separate(extra,
into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
sep = ",") %>%
select(gender:time) %>%
bind_rows(df.exp2.causal_first %>%
separate(extra,
into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
sep = ",") %>%
select(gender:time)) %>%
mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>%
mutate_if(is.character, ~ as.numeric(.)) %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) n age_mean age_sd n_female time_mean time_sd
1 82 35 12.2 45 21.25 5.11
# read model predictions
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "teleport")]
df.exp2.model = files %>%
map(~ read.csv(str_c(path, .))) %>%
reduce(rbind) %>%
set_names(c("clip", "noise", "prediction")) %>%
left_join(df.exp2.clipinfo, by = "clip") %>%
mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))
# calculate mean counterfactual judgments
df.exp2.counterfactual.means = df.exp2.combined.long %>%
filter(question == "counterfactual",
clip <= 18) %>%
# mutate(rating = abs(rating)) %>%
group_by(clip) %>%
summarize(rating = mean(rating)) %>%
ungroup() %>%
left_join(df.exp2.clipinfo, by = "clip")
# find noisy simulation model that correlates best with the mean counterfactual judgments
df.exp2.counterfactual.model = df.exp2.model %>%
group_by(noise) %>%
nest() %>%
mutate(correlation = map(data, ~ cor(.x$prediction,
df.exp2.counterfactual.means$rating)),
correlation = unlist(correlation)) %>%
ungroup() %>%
filter(correlation == max(correlation)) %>%
unnest(data) %>%
mutate(prediction = prediction * 100)# Model predictions based on counterfactual judgments
df.exp2.causal.model = df.exp2.combined.long %>%
filter(
question == "counterfactual",
clip <= 18
) %>%
group_by(clip, condition) %>%
summarize(rating = mean(rating)) %>%
spread(condition, rating) %>%
left_join(df.exp2.combined.long %>%
filter(
question == "counterfactual",
clip <= 18
) %>%
group_by(clip) %>%
summarize(combined = mean(rating)),
by = "clip") %>%
left_join(df.exp2.counterfactual.model,
by = "clip") %>%
select(-c(correlation, noise)) %>%
rename(simulation = prediction) %>%
mutate_at(.vars = vars(simulation, causal_first, counterfactual_first, combined), .funs = list(~ ifelse(outcome_actual == 1, 100 - ., .)))df.plot = df.exp2.combined.long %>%
filter(question == "counterfactual",
clip <= 18) %>%
mutate(clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
left_join(df.exp2.counterfactual.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
gather("model", "value", c(rating, prediction)) %>%
mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
colorindex = ifelse(model != "rating", 3, colorindex),
colorindex = as.factor(colorindex),
index_actual = factor(index_actual,
levels = c("actual miss",
"actual close",
"actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")),
model = factor(model,
levels = c("rating", "prediction")))
df.text = df.plot %>%
expand(index_actual, index_counterfactual, clipindex, model) %>%
mutate(
label = rep(c(
"1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
"7", "8 (C)", "9", "10", "11", "12 (C)",
"13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"
), each = 2),
label = ifelse(model != "rating", NA, label),
y = -15,
x = clipindex,
colorindex = NA
)
actual_counterfactual_plot(df.plot, "counterfactual judgment")
# ggsave("../../figures/plots/exp2_counterfactual_bars.pdf",
# width = 8,
# height = 6)df.exp2.counterfactual.means %>%
left_join(df.exp2.counterfactual.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
summarize(simulation_r = cor(rating, prediction),
simulation_rmse = rmse(rating, prediction)) %>%
mutate_all(~ round(., 2))# A tibble: 1 x 2
simulation_r simulation_rmse
<dbl> <dbl>
1 0.97 14.6
# random intercepts & slopes model
fit.brm_exp2_counterfactual = brm(formula = rating ~ 1 + condition * (index_actual + index_counterfactual) + (1 + index_actual + index_counterfactual | participant),
data = df.exp2.combined.long %>%
na.omit() %>%
filter(question == "counterfactual"),
cores = 2,
seed = 1,
file = "cache/brm_exp2_counterfactual")
fit.brm_exp2_counterfactual %>%
summary() Family: gaussian
Links: mu = identity; sigma = identity
Formula: rating ~ 1 + condition * (index_actual + index_counterfactual) + (1 + index_actual + index_counterfactual | participant)
Data: df.exp2.combined.long %>% na.omit() %>% filter(que (Number of observations: 1476)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~participant (Number of levels: 82)
Estimate
sd(Intercept) 8.80
sd(index_actualactualclose) 6.76
sd(index_actualactualhit) 2.37
sd(index_counterfactualcounterfactualclose) 14.65
sd(index_counterfactualcounterfactualhit) 22.12
cor(Intercept,index_actualactualclose) 0.13
cor(Intercept,index_actualactualhit) -0.02
cor(index_actualactualclose,index_actualactualhit) 0.18
cor(Intercept,index_counterfactualcounterfactualclose) -0.41
cor(index_actualactualclose,index_counterfactualcounterfactualclose) -0.04
cor(index_actualactualhit,index_counterfactualcounterfactualclose) -0.05
cor(Intercept,index_counterfactualcounterfactualhit) -0.86
cor(index_actualactualclose,index_counterfactualcounterfactualhit) -0.20
cor(index_actualactualhit,index_counterfactualcounterfactualhit) -0.11
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 0.36
Est.Error
sd(Intercept) 1.84
sd(index_actualactualclose) 3.07
sd(index_actualactualhit) 1.78
sd(index_counterfactualcounterfactualclose) 2.43
sd(index_counterfactualcounterfactualhit) 2.57
cor(Intercept,index_actualactualclose) 0.31
cor(Intercept,index_actualactualhit) 0.40
cor(index_actualactualclose,index_actualactualhit) 0.42
cor(Intercept,index_counterfactualcounterfactualclose) 0.20
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 0.29
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 0.38
cor(Intercept,index_counterfactualcounterfactualhit) 0.08
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 0.27
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 0.39
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 0.16
l-95% CI
sd(Intercept) 5.13
sd(index_actualactualclose) 0.58
sd(index_actualactualhit) 0.10
sd(index_counterfactualcounterfactualclose) 9.96
sd(index_counterfactualcounterfactualhit) 17.37
cor(Intercept,index_actualactualclose) -0.46
cor(Intercept,index_actualactualhit) -0.74
cor(index_actualactualclose,index_actualactualhit) -0.69
cor(Intercept,index_counterfactualcounterfactualclose) -0.73
cor(index_actualactualclose,index_counterfactualcounterfactualclose) -0.61
cor(index_actualactualhit,index_counterfactualcounterfactualclose) -0.73
cor(Intercept,index_counterfactualcounterfactualhit) -0.97
cor(index_actualactualclose,index_counterfactualcounterfactualhit) -0.70
cor(index_actualactualhit,index_counterfactualcounterfactualhit) -0.78
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 0.01
u-95% CI
sd(Intercept) 12.58
sd(index_actualactualclose) 12.43
sd(index_actualactualhit) 6.67
sd(index_counterfactualcounterfactualclose) 19.56
sd(index_counterfactualcounterfactualhit) 27.51
cor(Intercept,index_actualactualclose) 0.74
cor(Intercept,index_actualactualhit) 0.74
cor(index_actualactualclose,index_actualactualhit) 0.85
cor(Intercept,index_counterfactualcounterfactualclose) 0.02
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 0.55
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 0.71
cor(Intercept,index_counterfactualcounterfactualhit) -0.66
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 0.37
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 0.66
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 0.63
Rhat
sd(Intercept) 1.00
sd(index_actualactualclose) 1.01
sd(index_actualactualhit) 1.00
sd(index_counterfactualcounterfactualclose) 1.00
sd(index_counterfactualcounterfactualhit) 1.00
cor(Intercept,index_actualactualclose) 1.00
cor(Intercept,index_actualactualhit) 1.00
cor(index_actualactualclose,index_actualactualhit) 1.00
cor(Intercept,index_counterfactualcounterfactualclose) 1.00
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 1.00
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 1.01
cor(Intercept,index_counterfactualcounterfactualhit) 1.00
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 1.00
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 1.00
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 1.00
Bulk_ESS
sd(Intercept) 1720
sd(index_actualactualclose) 516
sd(index_actualactualhit) 1448
sd(index_counterfactualcounterfactualclose) 1135
sd(index_counterfactualcounterfactualhit) 1621
cor(Intercept,index_actualactualclose) 1418
cor(Intercept,index_actualactualhit) 3693
cor(index_actualactualclose,index_actualactualhit) 2473
cor(Intercept,index_counterfactualcounterfactualclose) 887
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 474
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 331
cor(Intercept,index_counterfactualcounterfactualhit) 652
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 1422
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 1323
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 2216
Tail_ESS
sd(Intercept) 1686
sd(index_actualactualclose) 639
sd(index_actualactualhit) 2198
sd(index_counterfactualcounterfactualclose) 1936
sd(index_counterfactualcounterfactualhit) 2274
cor(Intercept,index_actualactualclose) 2112
cor(Intercept,index_actualactualhit) 3102
cor(index_actualactualclose,index_actualactualhit) 2861
cor(Intercept,index_counterfactualcounterfactualclose) 1737
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 771
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 714
cor(Intercept,index_counterfactualcounterfactualhit) 1299
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 1430
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 2815
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 3157
Population-Level Effects:
Estimate
Intercept 13.38
conditioncounterfactual_first -1.11
index_actualactualclose -1.48
index_actualactualhit -4.42
index_counterfactualcounterfactualclose 37.94
index_counterfactualcounterfactualhit 73.80
conditioncounterfactual_first:index_actualactualclose -0.55
conditioncounterfactual_first:index_actualactualhit -0.85
conditioncounterfactual_first:index_counterfactualcounterfactualclose 1.23
conditioncounterfactual_first:index_counterfactualcounterfactualhit -0.02
Est.Error
Intercept 2.58
conditioncounterfactual_first 3.65
index_actualactualclose 2.78
index_actualactualhit 2.40
index_counterfactualcounterfactualclose 3.33
index_counterfactualcounterfactualhit 4.17
conditioncounterfactual_first:index_actualactualclose 3.84
conditioncounterfactual_first:index_actualactualhit 3.41
conditioncounterfactual_first:index_counterfactualcounterfactualclose 4.79
conditioncounterfactual_first:index_counterfactualcounterfactualhit 5.92
l-95% CI
Intercept 8.30
conditioncounterfactual_first -7.99
index_actualactualclose -6.79
index_actualactualhit -9.15
index_counterfactualcounterfactualclose 31.43
index_counterfactualcounterfactualhit 65.54
conditioncounterfactual_first:index_actualactualclose -8.29
conditioncounterfactual_first:index_actualactualhit -7.60
conditioncounterfactual_first:index_counterfactualcounterfactualclose -8.37
conditioncounterfactual_first:index_counterfactualcounterfactualhit -11.96
u-95% CI
Intercept 18.38
conditioncounterfactual_first 6.04
index_actualactualclose 3.95
index_actualactualhit 0.20
index_counterfactualcounterfactualclose 44.60
index_counterfactualcounterfactualhit 81.99
conditioncounterfactual_first:index_actualactualclose 6.79
conditioncounterfactual_first:index_actualactualhit 5.87
conditioncounterfactual_first:index_counterfactualcounterfactualclose 10.75
conditioncounterfactual_first:index_counterfactualcounterfactualhit 11.59
Rhat
Intercept 1.00
conditioncounterfactual_first 1.00
index_actualactualclose 1.00
index_actualactualhit 1.00
index_counterfactualcounterfactualclose 1.00
index_counterfactualcounterfactualhit 1.00
conditioncounterfactual_first:index_actualactualclose 1.00
conditioncounterfactual_first:index_actualactualhit 1.00
conditioncounterfactual_first:index_counterfactualcounterfactualclose 1.00
conditioncounterfactual_first:index_counterfactualcounterfactualhit 1.00
Bulk_ESS
Intercept 2268
conditioncounterfactual_first 1966
index_actualactualclose 2898
index_actualactualhit 3619
index_counterfactualcounterfactualclose 2118
index_counterfactualcounterfactualhit 1818
conditioncounterfactual_first:index_actualactualclose 2827
conditioncounterfactual_first:index_actualactualhit 3265
conditioncounterfactual_first:index_counterfactualcounterfactualclose 2006
conditioncounterfactual_first:index_counterfactualcounterfactualhit 1684
Tail_ESS
Intercept 2863
conditioncounterfactual_first 2564
index_actualactualclose 2790
index_actualactualhit 3303
index_counterfactualcounterfactualclose 2627
index_counterfactualcounterfactualhit 2490
conditioncounterfactual_first:index_actualactualclose 3028
conditioncounterfactual_first:index_actualactualhit 3250
conditioncounterfactual_first:index_counterfactualcounterfactualclose 2617
conditioncounterfactual_first:index_counterfactualcounterfactualhit 2692
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 26.57 0.56 25.50 27.72 1.00 2894 2705
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
fit.brm_exp2_counterfactual %>%
tidy() %>%
filter(str_detect(term, "b_")) %>%
mutate(term = str_remove(term, "b_"),
term = tolower(term)) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
rename(`lower 95% HDI` = lower,
`upper 95% HDI` = upper) %>%
xtable() %>%
print(include.rownames = F,
booktabs = T)% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:05:24 2019
\begin{table}[ht]
\centering
\begin{tabular}{lrrrr}
\toprule
term & estimate & std.error & lower 95\% HDI & upper 95\% HDI \\
\midrule
intercept & 13.38 & 2.58 & 9.12 & 17.55 \\
conditioncounterfactual\_first & -1.11 & 3.65 & -6.99 & 4.82 \\
index\_actualactualclose & -1.48 & 2.78 & -6.02 & 3.11 \\
index\_actualactualhit & -4.42 & 2.40 & -8.32 & -0.46 \\
index\_counterfactualcounterfactualclose & 37.94 & 3.33 & 32.39 & 43.35 \\
index\_counterfactualcounterfactualhit & 73.80 & 4.17 & 67.00 & 80.57 \\
conditioncounterfactual\_first:index\_actualactualclose & -0.55 & 3.84 & -7.11 & 5.73 \\
conditioncounterfactual\_first:index\_actualactualhit & -0.85 & 3.41 & -6.52 & 4.85 \\
conditioncounterfactual\_first:index\_counterfactualcounterfactualclose & 1.23 & 4.79 & -6.62 & 9.18 \\
conditioncounterfactual\_first:index\_counterfactualcounterfactualhit & -0.02 & 5.92 & -9.86 & 9.55 \\
\bottomrule
\end{tabular}
\end{table}
# which causal condition?
# condition.name = c("causal_first")
condition.name = c("counterfactual_first")
# condition.name = c("causal_first", "counterfactual_first")
# which model for the counterfactual simulations
# model.name = c("simulation")
# model.name = c("causal_first")
# model.name = c("counterfactual_first")
model.name = c("combined")
df.plot = df.exp2.combined.long %>%
filter(
question == "causal",
clip <= 18
) %>%
filter(condition %in% condition.name) %>%
mutate(
rating = abs(rating),
clipindex = rep(c(1, 2), nrow(.) / 2)
) %>%
left_join(df.exp2.causal.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
gather("model", "value", c("rating", model.name)) %>%
mutate(model = factor(model, levels = c("rating", model.name))) %>%
mutate(
colorindex = ifelse(outcome_actual == 1, 2, 1),
colorindex = ifelse(model != "rating", 3, colorindex),
colorindex = as.factor(colorindex),
index_actual = factor(index_actual,
levels = c("actual miss",
"actual close",
"actual hit")),
index_counterfactual = factor(index_counterfactual,
levels = c("counterfactual miss",
"counterfactual close",
"counterfactual hit")))
df.text = df.plot %>%
expand(index_actual, index_counterfactual, clipindex, model) %>%
mutate(
label = rep(c(
"1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
"7", "8 (C)", "9", "10", "11", "12 (C)",
"13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"
), each = 2),
label = ifelse(model != "rating", NA, label),
y = -15,
x = clipindex,
colorindex = NA
)
actual_counterfactual_plot(df.plot, "causal judgment")Warning: Removed 738 rows containing missing values (geom_point).
# ggsave(str_c("../../figures/plots/exp2_", condition.name ,"_causal_bars.pdf"),
# width = 8,
# height = 6)df.exp2.combined.long %>%
filter(question == "causal",
clip <= 18) %>%
mutate(rating = abs(rating)) %>%
group_by(clip,
outcome_actual,
outcome_counterfactual,
index_actual,
index_counterfactual) %>%
summarize(rating = mean(rating)) %>%
ungroup() %>%
left_join(df.exp2.causal.model,
by = c("clip",
"outcome_actual",
"outcome_counterfactual",
"index_actual",
"index_counterfactual")) %>%
summarize(empirical_r = cor(rating, combined),
empirical_rmse = rmse(rating, combined),
simulation_r = cor(rating, simulation),
simulation_rmse = rmse(rating, simulation)) %>%
mutate_all(~ round(., 2))# A tibble: 1 x 4
empirical_r empirical_rmse simulation_r simulation_rmse
<dbl> <dbl> <dbl> <dbl>
1 0.97 8.36 0.94 19.4
# random intercepts + slopes model
fit.brm_exp2_causal = brm(formula = rating ~ 1 + condition * (index_actual + index_counterfactual + outcome_actual) +
(1 + index_actual + index_counterfactual + outcome_actual | participant),
data = df.exp2.combined.long %>%
filter(question == "causal",
clip <= 18),
cores = 2,
seed = 1,
file = "cache/brm_exp2_causal")
fit.brm_exp2_causal %>%
summary() Family: gaussian
Links: mu = identity; sigma = identity
Formula: rating ~ 1 + condition * (index_actual + index_counterfactual + outcome_actual) + (1 + index_actual + index_counterfactual + outcome_actual | participant)
Data: df.exp2.combined.long %>% filter(question == "caus (Number of observations: 1476)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~participant (Number of levels: 82)
Estimate
sd(Intercept) 15.12
sd(index_actualactualclose) 3.57
sd(index_actualactualhit) 6.72
sd(index_counterfactualcounterfactualclose) 15.91
sd(index_counterfactualcounterfactualhit) 22.85
sd(outcome_actual) 23.85
cor(Intercept,index_actualactualclose) -0.16
cor(Intercept,index_actualactualhit) -0.43
cor(index_actualactualclose,index_actualactualhit) 0.18
cor(Intercept,index_counterfactualcounterfactualclose) 0.02
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 0.19
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 0.03
cor(Intercept,index_counterfactualcounterfactualhit) -0.62
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 0.18
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 0.39
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 0.20
cor(Intercept,outcome_actual) -0.69
cor(index_actualactualclose,outcome_actual) 0.09
cor(index_actualactualhit,outcome_actual) 0.02
cor(index_counterfactualcounterfactualclose,outcome_actual) -0.13
cor(index_counterfactualcounterfactualhit,outcome_actual) 0.07
Est.Error
sd(Intercept) 2.20
sd(index_actualactualclose) 2.41
sd(index_actualactualhit) 3.79
sd(index_counterfactualcounterfactualclose) 2.75
sd(index_counterfactualcounterfactualhit) 2.88
sd(outcome_actual) 2.96
cor(Intercept,index_actualactualclose) 0.35
cor(Intercept,index_actualactualhit) 0.31
cor(index_actualactualclose,index_actualactualhit) 0.37
cor(Intercept,index_counterfactualcounterfactualclose) 0.19
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 0.34
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 0.32
cor(Intercept,index_counterfactualcounterfactualhit) 0.11
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 0.35
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 0.31
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 0.18
cor(Intercept,outcome_actual) 0.12
cor(index_actualactualclose,outcome_actual) 0.34
cor(index_actualactualhit,outcome_actual) 0.30
cor(index_counterfactualcounterfactualclose,outcome_actual) 0.18
cor(index_counterfactualcounterfactualhit,outcome_actual) 0.17
l-95% CI
sd(Intercept) 11.02
sd(index_actualactualclose) 0.16
sd(index_actualactualhit) 0.48
sd(index_counterfactualcounterfactualclose) 10.82
sd(index_counterfactualcounterfactualhit) 17.57
sd(outcome_actual) 18.35
cor(Intercept,index_actualactualclose) -0.75
cor(Intercept,index_actualactualhit) -0.87
cor(index_actualactualclose,index_actualactualhit) -0.59
cor(Intercept,index_counterfactualcounterfactualclose) -0.34
cor(index_actualactualclose,index_counterfactualcounterfactualclose) -0.53
cor(index_actualactualhit,index_counterfactualcounterfactualclose) -0.61
cor(Intercept,index_counterfactualcounterfactualhit) -0.80
cor(index_actualactualclose,index_counterfactualcounterfactualhit) -0.56
cor(index_actualactualhit,index_counterfactualcounterfactualhit) -0.36
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) -0.17
cor(Intercept,outcome_actual) -0.88
cor(index_actualactualclose,outcome_actual) -0.60
cor(index_actualactualhit,outcome_actual) -0.56
cor(index_counterfactualcounterfactualclose,outcome_actual) -0.47
cor(index_counterfactualcounterfactualhit,outcome_actual) -0.26
u-95% CI
sd(Intercept) 19.73
sd(index_actualactualclose) 8.90
sd(index_actualactualhit) 14.42
sd(index_counterfactualcounterfactualclose) 21.60
sd(index_counterfactualcounterfactualhit) 28.81
sd(outcome_actual) 30.03
cor(Intercept,index_actualactualclose) 0.56
cor(Intercept,index_actualactualhit) 0.32
cor(index_actualactualclose,index_actualactualhit) 0.80
cor(Intercept,index_counterfactualcounterfactualclose) 0.41
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 0.78
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 0.64
cor(Intercept,index_counterfactualcounterfactualhit) -0.37
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 0.78
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 0.86
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 0.52
cor(Intercept,outcome_actual) -0.42
cor(index_actualactualclose,outcome_actual) 0.71
cor(index_actualactualhit,outcome_actual) 0.61
cor(index_counterfactualcounterfactualclose,outcome_actual) 0.23
cor(index_counterfactualcounterfactualhit,outcome_actual) 0.38
Rhat
sd(Intercept) 1.00
sd(index_actualactualclose) 1.00
sd(index_actualactualhit) 1.00
sd(index_counterfactualcounterfactualclose) 1.00
sd(index_counterfactualcounterfactualhit) 1.00
sd(outcome_actual) 1.00
cor(Intercept,index_actualactualclose) 1.00
cor(Intercept,index_actualactualhit) 1.00
cor(index_actualactualclose,index_actualactualhit) 1.00
cor(Intercept,index_counterfactualcounterfactualclose) 1.00
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 1.01
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 1.01
cor(Intercept,index_counterfactualcounterfactualhit) 1.00
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 1.01
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 1.02
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 1.00
cor(Intercept,outcome_actual) 1.00
cor(index_actualactualclose,outcome_actual) 1.00
cor(index_actualactualhit,outcome_actual) 1.01
cor(index_counterfactualcounterfactualclose,outcome_actual) 1.00
cor(index_counterfactualcounterfactualhit,outcome_actual) 1.01
Bulk_ESS
sd(Intercept) 1772
sd(index_actualactualclose) 910
sd(index_actualactualhit) 1024
sd(index_counterfactualcounterfactualclose) 1102
sd(index_counterfactualcounterfactualhit) 1360
sd(outcome_actual) 1221
cor(Intercept,index_actualactualclose) 3259
cor(Intercept,index_actualactualhit) 2183
cor(index_actualactualclose,index_actualactualhit) 1812
cor(Intercept,index_counterfactualcounterfactualclose) 1238
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 311
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 507
cor(Intercept,index_counterfactualcounterfactualhit) 1487
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 330
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 363
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 1388
cor(Intercept,outcome_actual) 881
cor(index_actualactualclose,outcome_actual) 556
cor(index_actualactualhit,outcome_actual) 689
cor(index_counterfactualcounterfactualclose,outcome_actual) 1250
cor(index_counterfactualcounterfactualhit,outcome_actual) 1664
Tail_ESS
sd(Intercept) 2323
sd(index_actualactualclose) 1374
sd(index_actualactualhit) 1620
sd(index_counterfactualcounterfactualclose) 1926
sd(index_counterfactualcounterfactualhit) 2216
sd(outcome_actual) 2105
cor(Intercept,index_actualactualclose) 2034
cor(Intercept,index_actualactualhit) 2671
cor(index_actualactualclose,index_actualactualhit) 2987
cor(Intercept,index_counterfactualcounterfactualclose) 2292
cor(index_actualactualclose,index_counterfactualcounterfactualclose) 590
cor(index_actualactualhit,index_counterfactualcounterfactualclose) 1166
cor(Intercept,index_counterfactualcounterfactualhit) 2044
cor(index_actualactualclose,index_counterfactualcounterfactualhit) 784
cor(index_actualactualhit,index_counterfactualcounterfactualhit) 872
cor(index_counterfactualcounterfactualclose,index_counterfactualcounterfactualhit) 2445
cor(Intercept,outcome_actual) 1398
cor(index_actualactualclose,outcome_actual) 1563
cor(index_actualactualhit,outcome_actual) 1197
cor(index_counterfactualcounterfactualclose,outcome_actual) 2261
cor(index_counterfactualcounterfactualhit,outcome_actual) 2723
Population-Level Effects:
Estimate
Intercept -24.07
conditioncounterfactual_first 11.47
index_actualactualclose 8.17
index_actualactualhit 13.42
index_counterfactualcounterfactualclose -30.20
index_counterfactualcounterfactualhit -50.29
outcome_actual 98.53
conditioncounterfactual_first:index_actualactualclose -5.39
conditioncounterfactual_first:index_actualactualhit -16.97
conditioncounterfactual_first:index_counterfactualcounterfactualclose -8.13
conditioncounterfactual_first:index_counterfactualcounterfactualhit -25.51
conditioncounterfactual_first:outcome_actual 10.72
Est.Error
Intercept 3.50
conditioncounterfactual_first 4.91
index_actualactualclose 3.44
index_actualactualhit 4.90
index_counterfactualcounterfactualclose 3.72
index_counterfactualcounterfactualhit 4.60
outcome_actual 5.40
conditioncounterfactual_first:index_actualactualclose 4.94
conditioncounterfactual_first:index_actualactualhit 7.02
conditioncounterfactual_first:index_counterfactualcounterfactualclose 5.17
conditioncounterfactual_first:index_counterfactualcounterfactualhit 6.44
conditioncounterfactual_first:outcome_actual 7.60
l-95% CI
Intercept -30.87
conditioncounterfactual_first 1.62
index_actualactualclose 1.47
index_actualactualhit 3.74
index_counterfactualcounterfactualclose -37.61
index_counterfactualcounterfactualhit -59.57
outcome_actual 87.81
conditioncounterfactual_first:index_actualactualclose -15.00
conditioncounterfactual_first:index_actualactualhit -30.40
conditioncounterfactual_first:index_counterfactualcounterfactualclose -18.18
conditioncounterfactual_first:index_counterfactualcounterfactualhit -38.16
conditioncounterfactual_first:outcome_actual -4.16
u-95% CI
Intercept -17.24
conditioncounterfactual_first 21.22
index_actualactualclose 14.79
index_actualactualhit 22.69
index_counterfactualcounterfactualclose -23.19
index_counterfactualcounterfactualhit -41.47
outcome_actual 108.98
conditioncounterfactual_first:index_actualactualclose 4.35
conditioncounterfactual_first:index_actualactualhit -3.33
conditioncounterfactual_first:index_counterfactualcounterfactualclose 2.17
conditioncounterfactual_first:index_counterfactualcounterfactualhit -12.98
conditioncounterfactual_first:outcome_actual 26.02
Rhat
Intercept 1.00
conditioncounterfactual_first 1.00
index_actualactualclose 1.00
index_actualactualhit 1.00
index_counterfactualcounterfactualclose 1.00
index_counterfactualcounterfactualhit 1.00
outcome_actual 1.00
conditioncounterfactual_first:index_actualactualclose 1.00
conditioncounterfactual_first:index_actualactualhit 1.00
conditioncounterfactual_first:index_counterfactualcounterfactualclose 1.00
conditioncounterfactual_first:index_counterfactualcounterfactualhit 1.00
conditioncounterfactual_first:outcome_actual 1.00
Bulk_ESS
Intercept 1592
conditioncounterfactual_first 1214
index_actualactualclose 2313
index_actualactualhit 1993
index_counterfactualcounterfactualclose 2457
index_counterfactualcounterfactualhit 1597
outcome_actual 1983
conditioncounterfactual_first:index_actualactualclose 2251
conditioncounterfactual_first:index_actualactualhit 1990
conditioncounterfactual_first:index_counterfactualcounterfactualclose 2373
conditioncounterfactual_first:index_counterfactualcounterfactualhit 1212
conditioncounterfactual_first:outcome_actual 1712
Tail_ESS
Intercept 2398
conditioncounterfactual_first 1554
index_actualactualclose 2743
index_actualactualhit 2073
index_counterfactualcounterfactualclose 2375
index_counterfactualcounterfactualhit 1927
outcome_actual 2818
conditioncounterfactual_first:index_actualactualclose 2511
conditioncounterfactual_first:index_actualactualhit 2024
conditioncounterfactual_first:index_counterfactualcounterfactualclose 2225
conditioncounterfactual_first:index_counterfactualcounterfactualhit 2082
conditioncounterfactual_first:outcome_actual 2341
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 29.90 0.65 28.65 31.18 1.00 2704 3029
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
df.exp2.causal.posterior = df.exp2.combined.long %>%
filter(question == "causal",
clip <= 18) %>%
group_by(clip, index_actual, index_counterfactual, outcome_actual, condition) %>%
summarize(value = mean(rating)) %>%
ungroup() %>%
add_fitted_draws(newdata = .,
model = fit.brm_exp2_causal,
re_formula = NA) %>%
ungroup() %>%
select(clip, condition, .value, .draw) %>%
spread(clip, .value)
# difference between clips 13 and 14 versus 8
df.exp2.causal.posterior %>%
mutate(difference = (`13` + `14`) / 2 - `8`) %>%
pull(difference) %>%
mean_hdi() %>%
mutate_if(is.numeric, ~round(., 2)) %>%
summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])"))
# difference between clips 5 and 6 versus 11
df.exp2.causal.posterior %>%
mutate(difference = (`5` + `6`) / 2 - `11`) %>%
pull(difference) %>%
mean_hdi() %>%
mutate_if(is.numeric, ~round(., 2)) %>%
summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])")) difference
1 (-0.55 [-12.09, 10.97])
difference
1 (-5.48 [-13.8, 2.73])
fit.brm_exp2_causal %>%
tidy() %>%
filter(str_detect(term, "b_")) %>%
mutate(term = str_remove(term, "b_"),
term = tolower(term)) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
rename(`lower 95% HDI` = lower,
`upper 95% HDI` = upper) %>%
xtable() %>%
print(include.rownames = F,
booktabs = T)% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:05:29 2019
\begin{table}[ht]
\centering
\begin{tabular}{lrrrr}
\toprule
term & estimate & std.error & lower 95\% HDI & upper 95\% HDI \\
\midrule
intercept & -24.07 & 3.50 & -29.80 & -18.27 \\
conditioncounterfactual\_first & 11.47 & 4.91 & 3.29 & 19.40 \\
index\_actualactualclose & 8.17 & 3.44 & 2.39 & 13.62 \\
index\_actualactualhit & 13.42 & 4.90 & 5.42 & 21.30 \\
index\_counterfactualcounterfactualclose & -30.20 & 3.72 & -36.41 & -24.36 \\
index\_counterfactualcounterfactualhit & -50.29 & 4.60 & -57.94 & -42.84 \\
outcome\_actual & 98.53 & 5.40 & 89.53 & 107.57 \\
conditioncounterfactual\_first:index\_actualactualclose & -5.39 & 4.94 & -13.38 & 2.80 \\
conditioncounterfactual\_first:index\_actualactualhit & -16.97 & 7.02 & -28.61 & -5.75 \\
conditioncounterfactual\_first:index\_counterfactualcounterfactualclose & -8.13 & 5.17 & -16.51 & 0.35 \\
conditioncounterfactual\_first:index\_counterfactualcounterfactualhit & -25.51 & 6.44 & -36.53 & -15.18 \\
conditioncounterfactual\_first:outcome\_actual & 10.72 & 7.60 & -1.63 & 23.34 \\
\bottomrule
\end{tabular}
\end{table}
# clipinfo
df.exp3.clipinfo = tibble(
clip = 1:32,
outcome_both = c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1),
outcome_a = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1),
outcome_b = c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1),
outcome_none = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
clipindex = rep(1:2, 16))
# COUNTERFACTUAL JUDGMENTS
df.exp3.counterfactual = read.csv("../../data/experiment3_counterfactual.csv", header = T, stringsAsFactors = F)
# demographics
df.exp3.counterfactual.demographics = df.exp3.counterfactual$extra %>%
strsplit("\n") %>%
unlist() %>%
strsplit(",") %>%
unlist() %>%
matrix(ncol = 6, byrow = T) %>%
as.data.frame(stringsAsFactors = F) %>%
mutate_all("as.numeric") %>%
setNames(c("gender", "age", "difficulty", "smoothness", "time", "condition")) %>%
mutate(
participant = 1:nrow(.),
condition = factor(condition, levels = 20:21, labels = c("A", "B")),
gender = factor(gender, levels = 1:2, labels = c("female", "male"))
) %>%
select(participant, condition, gender, age, time, smoothness, difficulty)
# judgments
df.exp3.counterfactual.long = df.exp3.counterfactual$data %>%
strsplit("\n") %>%
unlist() %>%
strsplit(",") %>%
unlist() %>%
matrix(ncol = 64, byrow = T) %>%
as.data.frame(stringsAsFactors = F) %>%
mutate_all("as.numeric") %>%
setNames(str_c(c("clip_", "rating_"), rep(1:32, each = 2))) %>%
mutate(participant = 1:nrow(.)) %>%
gather("index", "value", -participant) %>%
separate(index, into = c("index", "order")) %>%
mutate(order = as.numeric(order)) %>%
spread(index, value) %>%
left_join(select(df.exp3.counterfactual.demographics, participant, ball = condition),
by = "participant") %>%
select(participant, ball, clip, everything()) %>%
arrange(participant, clip)
# CAUSAL JUDGMENTS
df.exp3.causal = read.csv("../../data/experiment3_causal.csv", header = T, stringsAsFactors = F)
# demographics
df.exp3.causal.demographics = df.exp3.causal$extra %>%
strsplit("\n") %>%
unlist() %>%
strsplit(",") %>%
unlist() %>%
matrix(nrow = nrow(df.exp3.causal), byrow = T) %>%
as.data.frame(stringsAsFactors = F) %>%
mutate_all("as.numeric") %>%
setNames(c("gender", "age", "difficulty", "smoothness", "time", "counterbalance", "replayType", "experiment")) %>%
mutate(participant = 1:nrow(.),
counterbalance = ifelse(participant %in% c(1, 3, 5, 6, 9), 1, counterbalance),
gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>%
select(participant, gender, age, counterbalance, time, smoothness, difficulty)
# judgments
df.exp3.causal.long = df.exp3.causal$data %>%
strsplit("\n") %>%
unlist() %>%
strsplit(",") %>%
unlist() %>%
matrix(ncol = 192, byrow = T) %>%
as.data.frame(stringsAsFactors = F) %>%
mutate_all("as.numeric") %>%
setNames(str_c(c("clip_", "X_", "Y_", "Z_", "rating1_", "rating2_"), rep(1:32, each = 6))) %>%
select(contains("clip"), contains("rating")) %>%
mutate(participant = 1:nrow(.)) %>%
gather("index", "value", -participant) %>%
separate(index, into = c("index", "order")) %>%
mutate(order = as.numeric(order)) %>%
spread(index, value) %>%
left_join(select(df.exp3.causal.demographics, participant, counterbalance),
by = "participant") %>%
arrange(participant, clip) %>%
mutate(
A = ifelse(counterbalance == 1, rating1, rating2),
B = ifelse(counterbalance == 1, rating2, rating1)
) %>%
select(-rating1, -rating2) %>%
gather("ball", "rating", c(A, B)) %>%
select(participant, clip, ball, rating, order) %>%
arrange(participant, clip, ball)# counterfactual condition
df.exp3.counterfactual.demographics %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2))
# causal condition
df.exp3.causal.demographics %>%
summarize(n = nrow(.),
age_mean = mean(age) %>% round(0),
age_sd = sd(age) %>% round(1),
n_female = sum(gender == "female"),
time_mean = mean(time) %>% round(2),
time_sd = sd(time) %>% round(2)) n age_mean age_sd n_female time_mean time_sd
1 80 33 10.1 34 18.08 4.63
n age_mean age_sd n_female time_mean time_sd
1 41 34 10.5 21 21.19 4.96
# read model predictions
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "3ball")]
df.exp3.model = files %>%
map(~ read.csv(str_c(path, .))) %>%
reduce(rbind) %>%
rename(clip = trial)# calculate mean counterfactual judgments
df.exp3.counterfactual.means = df.exp3.counterfactual.long %>%
group_by(clip, ball) %>%
summarize(rating_mean = mean(rating),
rating_low = smean.cl.boot(rating)[2],
rating_high = smean.cl.boot(rating)[3]) %>%
ungroup()
# find noisy simulation model that correlates best with the mean counterfactual judgments
df.exp3.model.counterfactual = df.exp3.model %>%
select(clip, contains("whether"), noise) %>%
gather("ball", "prediction", c(A_whether, B_whether)) %>%
mutate(ball = str_remove(ball, "_whether")) %>%
arrange(noise, clip, ball) %>%
left_join(df.exp3.clipinfo, by = "clip") %>%
mutate(prediction = ifelse(outcome_both == 1, 1 - prediction, prediction)) %>%
group_by(noise) %>%
nest() %>%
mutate(correlation = map(data, ~ cor(.x$prediction,
df.exp3.counterfactual.means$rating_mean)),
correlation = unlist(correlation)) %>%
ungroup() %>%
filter(correlation == max(correlation)) %>%
unnest(data) %>%
mutate(prediction = prediction * 100)# calculate mean causal judgments
df.exp3.causal.means = df.exp3.causal.long %>%
group_by(clip, ball) %>%
summarize(rating_mean = mean(rating),
rating_low = smean.cl.boot(rating)[2],
rating_high = smean.cl.boot(rating)[3]) %>%
ungroup()
df.exp3.model.causal = df.exp3.model %>%
# take into account difference-making
mutate_at(.vars = vars(A_how:A_robust),
.funs = list(~ . * A_difference)) %>%
mutate_at(.vars = vars(B_how:B_robust),
.funs = list(~ . * B_difference)) %>%
# choose model based on best fit with counterfactual data
filter(noise == unique(df.exp3.model.counterfactual$noise)) %>%
gather("index", "value", A_difference:B_robust) %>%
separate(index, into = c("ball", "aspect")) %>%
spread(aspect, value) %>%
select(clip, ball:whether)l.features = fromJSON("data/features.json")
df.features.initial = l.features[["appearance"]] %>%
cbind(l.features[["initialVelocity"]][["ballA"]]) %>%
cbind(l.features[["initialVelocity"]][["ballB"]]) %>%
cbind(l.features[["initialVelocity"]][["ballE"]]) %>%
setNames(c(
paste0("ball", LETTERS[c(1, 2, 5)], "_appearance"),
"ballA_velx",
"ballA_vely",
"ballB_velx",
"ballB_vely",
"ballE_velx",
"ballE_vely"
)) %>%
mutate(clip = 1:nrow(.)) %>%
pivot_longer(cols = starts_with("ball")) %>%
separate(name, into = c("ball", "property")) %>%
pivot_wider(names_from = property, values_from = value) %>%
mutate(ball = str_remove_all(ball, pattern = "ball"))
# information about collisions
df.features.collisions = data.frame()
ballnames = c("ballA", "ballB", "ballE")
for (i in 1:32) {
df.tmp = data.frame()
for (j in 1:length(ballnames)) {
ncollisions = length(l.features[["collisions"]][[ballnames[j]]][["object"]][[i]])
if (ncollisions > 0) {
tmp = data.frame(
ball = ballnames[j],
object = l.features[["collisions"]][[ballnames[j]]][["object"]][[i]],
time = l.features[["collisions"]][[ballnames[j]]][["time"]][[i]],
pre.velx = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["x"]],
pre.vely = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["y"]],
post.velx = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["x"]],
post.vely = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["y"]],
ncollision = 1:ncollisions
)
df.tmp = df.tmp %>%
rbind(tmp)
}
}
df.features.collisions = df.features.collisions %>%
rbind(df.tmp %>%
mutate(clip = i) %>%
select(clip, ball, everything()))
}
# find end of clip
tmp = df.features.collisions %>%
group_by(clip) %>%
filter(
ball == "ballE",
str_detect(object, "goal|Left|Right")
) %>%
group_by(clip) %>%
mutate(endclip = max(time)) %>%
select(clip, endclip) %>%
ungroup() %>%
distinct()
# remove events that happen after the end of the clip
df.features.collisions = df.features.collisions %>%
left_join(tmp) %>%
mutate(endclip = ifelse(is.na(endclip), 400, endclip)) %>%
group_by(clip) %>%
filter(time <= endclip)Joining, by = "clip"
# E initially at rest or moving?
tmp.movingE = df.features.initial %>%
filter(ball == "E") %>%
mutate(E.moving = ifelse(velx == 0 & vely == 0, 1, 0)) %>%
select(clip, E.moving)
# contact with ball E
tmp.contactE = df.features.collisions %>%
filter(object == "ballE") %>%
mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
mutate(contactE = 1) %>%
select(clip, ball, contactE) %>%
group_by(clip, ball) %>%
summarise(contactE = any(contactE %>% as.logical()) * 1) %>%
left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .) %>%
mutate(contactE = ifelse(is.na(contactE), 0, contactE))Joining, by = c("clip", "ball")
# number of collisions
tmp.ncollisions = df.features.collisions %>%
filter(
str_detect(object, "ball"),
ball != "ballE"
) %>%
mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
group_by(clip, ball) %>%
count() %>%
rename(ncollision = n)
# exclusivity (no contact with ball E by the other ball)
tmp.exclusive = df.features.collisions %>%
filter(
str_detect(object, "ball"),
ball != "ballE"
) %>%
count(clip, ball, object) %>%
mutate(
AE = ifelse(ball == "ballA" & object == "ballE", 1, 0),
BE = ifelse(ball == "ballB" & object == "ballE", 1, 0)
) %>%
group_by(clip) %>%
summarise(
AE = any(AE %>% as.logical()) * 1,
BE = any(BE %>% as.logical()) * 1
) %>%
mutate(
A = ifelse(AE == 1 & BE == 0, 1, 0),
B = ifelse(AE == 0 & BE == 1, 1, 0)
) %>%
select(clip, A, B) %>%
gather(ball, exclusive, -clip)
# impact
func_angle = function(x, y) {
dot.prod = x %*% y
norm.x = norm(x, type = "2")
norm.y = norm(y, type = "2")
theta = acos(dot.prod / (norm.x * norm.y))
as.numeric(theta)
}
# impact on ball E
tmp.impactE = df.features.collisions %>%
rowwise() %>%
filter(
str_detect(object, "ball"),
ball == "ballE"
) %>%
mutate(
pre.speed = abs(pre.velx) + abs(pre.vely),
post.speed = abs(post.velx) + abs(post.vely),
speed.diff = post.speed - pre.speed,
direction.diff = func_angle(c(pre.velx, pre.vely), c(post.velx, post.vely)),
direction.diff = ifelse(is.na(direction.diff), abs(atan2(post.vely - pre.vely, post.velx - pre.velx)), direction.diff)
) %>%
ungroup() %>%
group_by(clip, object) %>%
summarise(
speed.diff = sum(speed.diff),
direction.diff = sum(direction.diff)
) %>%
rename(ball = object) %>%
mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .) %>%
mutate_at(.vars = vars(speed.diff, direction.diff), .funs = ~ ifelse(is.na(.), 0, .)) %>%
arrange(clip, ball) %>%
rename(
E.speed.diff = speed.diff,
E.direction.diff = direction.diff
)Joining, by = c("clip", "ball")
# total impact
tmp.impactTotal = df.features.collisions %>%
rowwise() %>%
filter(str_detect(object, "ballA|ballB")) %>%
mutate(
pre.speed = abs(pre.velx) + abs(pre.vely),
post.speed = abs(post.velx) + abs(post.vely),
speed.diff = post.speed - pre.speed,
direction.diff = func_angle(c(pre.velx, pre.vely), c(post.velx, post.vely)),
direction.diff = ifelse(is.na(direction.diff), abs(atan2(post.vely - pre.vely, post.velx - pre.velx)), direction.diff)
) %>%
ungroup() %>%
group_by(clip, object) %>%
summarise(
speed.diff = sum(speed.diff),
direction.diff = sum(direction.diff)
) %>%
rename(ball = object) %>%
mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .) %>%
mutate_at(.vars = vars(speed.diff, direction.diff), .funs = ~ ifelse(is.na(.), 0, .)) %>%
arrange(clip, ball) %>%
rename(
total.speed.diff = speed.diff,
total.direction.diff = direction.diff
)Joining, by = c("clip", "ball")
# transfer of force
tmp.transfer = df.features.collisions %>%
filter(
str_detect(ball, "ballA|ballB"),
str_detect(object, "ball")
) %>%
mutate(
AE = ifelse(ball == "ballA" & object == "ballE", 1, NA),
BE = ifelse(ball == "ballB" & object == "ballE", 1, NA),
AB = ifelse((ball == "ballA" & object == "ballB"), 1, NA)
) %>%
mutate_at(.vars = vars(AE, BE, AB), .funs = ~ . * time) %>%
filter(!is.na(AE) | !is.na(BE) | !is.na(AB)) %>%
arrange(clip, time) %>%
group_by(clip) %>%
summarise(
A = any(!is.na(AE)),
B = any(!is.na(BE)),
A = ifelse(any(!is.na(AB)) & any(!is.na(BE)) & max(AB, na.rm = T) < min(BE, na.rm = T), T, A),
B = ifelse(any(!is.na(AB)) & any(!is.na(AE)) & max(AB, na.rm = T) < min(AE, na.rm = T), T, B)
) %>%
gather(ball, transfer, -clip) %>%
arrange(clip, ball)Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in min(BE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf
Warning in max(AB, na.rm = T): no non-missing arguments to max; returning -
Inf
Warning in min(AE, na.rm = T): no non-missing arguments to min; returning
Inf
# collect features
df.features = df.features.initial %>%
filter(ball != "E") %>%
mutate(ball = factor(ball)) %>%
mutate(
moving = ifelse(velx == 0 & vely == 0, 0, 1),
speed = abs(velx) + abs(vely)
) %>%
left_join(tmp.contactE) %>%
left_join(tmp.impactE) %>%
left_join(tmp.impactTotal) %>%
left_join(tmp.transfer %>% mutate(transfer = transfer * 1)) %>%
left_join(tmp.movingE) %>%
left_join(tmp.exclusive) %>%
select(-c(appearance, velx, vely))Joining, by = c("clip", "ball")
Joining, by = c("clip", "ball")
Joining, by = c("clip", "ball")
Joining, by = c("clip", "ball")
Warning: Column `ball` joining factor and character vector, coercing into
character vector
Joining, by = "clip"
Joining, by = c("clip", "ball")
df.plot = df.exp3.counterfactual.long %>%
left_join(df.exp3.model.counterfactual %>%
select(clip, ball, prediction) %>%
mutate(ball = as.factor(ball)),
by = c("clip", "ball")
) %>%
left_join(df.exp3.clipinfo,
by = "clip") %>%
gather("model", "value", c(rating, prediction)) %>%
mutate(
colorindex = 1,
colorindex = ifelse(model == "rating" & ball == "A" & outcome_b == 1, 2, colorindex),
colorindex = ifelse(model == "rating" & ball == "B" & outcome_a == 1, 2, colorindex),
colorindex = ifelse(model == "prediction", 3, colorindex),
colorindex = as.factor(colorindex),
model = factor(model, levels = c("rating", "prediction"))
) %>%
arrange(participant, clip, ball)
df.text = df.plot %>%
expand(ball, clip, model) %>%
mutate(
label = ifelse(model == "rating" & ball == "A", clip, NA),
y = 110,
x = 1.7,
colorindex = NA
)
actual_counterfactual_threeballs_plot(df.plot, "counterfactual judgment")
# ggsave("../../figures/plots/exp3_counterfactual_bars.pdf",
# width = 12,
# height = 6)# best fitting model
df.exp3.counterfactual.means %>%
left_join(df.exp3.model.counterfactual) %>%
summarize(simulation_r = cor(rating_mean, prediction),
simulation_rmse = rmse(rating_mean, prediction)) %>%
mutate_all(~ round(., 2))Joining, by = c("clip", "ball")
Warning: Column `ball` joining factor and character vector, coercing into
character vector
# deterministic model
df.exp3.counterfactual.means %>%
left_join(df.exp3.clipinfo %>%
select(clip, outcome_a, outcome_b) %>%
gather("ball", "outcome", -clip) %>%
mutate(ball = factor(ball,
levels = c("outcome_a", "outcome_b"),
labels = c("B", "A")))) %>%
mutate(outcome = outcome * 100) %>%
summarize(simulation_r = cor(rating_mean, outcome),
simulation_rmse = rmse(rating_mean, outcome)) %>%
mutate_all(~ round(., 2))Joining, by = c("clip", "ball")
Warning: Column `ball` joining factors with different levels, coercing to
character vector
# A tibble: 1 x 2
simulation_r simulation_rmse
<dbl> <dbl>
1 0.87 20.7
# A tibble: 1 x 2
simulation_r simulation_rmse
<dbl> <dbl>
1 0.83 30.3
df.data = df.exp3.causal.long %>%
left_join(df.exp3.model.causal,
by = c("clip", "ball"))
fit.brm_exp3_causal_w = brm(formula = rating ~ 1 + whether + (1 + whether | participant),
data = df.data,
save_all_pars = T,
cores = 2,
seed = 1,
file = "cache/brm_exp3_causal_w")
fit.brm_exp3_causal_wh = brm(formula = rating ~ 1 + whether + how + (1 + whether + how | participant),
data = df.data,
save_all_pars = T,
cores = 2,
seed = 1,
file = "cache/brm_exp3_causal_wh")
fit.brm_exp3_causal_whs = brm(formula = rating ~ 1 + whether + how + sufficient + (1 + whether + how + sufficient | participant),
data = df.data,
save_all_pars = T,
cores = 2,
seed = 1,
file = "cache/brm_exp3_causal_whs")
fit.brm_exp3_causal_w %>% summary()
fit.brm_exp3_causal_wh %>% summary()
fit.brm_exp3_causal_whs %>% summary() Family: gaussian
Links: mu = identity; sigma = identity
Formula: rating ~ 1 + whether + (1 + whether | participant)
Data: df.data (Number of observations: 2624)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~participant (Number of levels: 41)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(Intercept) 7.76 1.64 4.78 11.03 1.00 1762
sd(whether) 8.39 2.93 2.47 13.90 1.01 677
cor(Intercept,whether) 0.00 0.37 -0.59 0.81 1.00 1100
Tail_ESS
sd(Intercept) 2343
sd(whether) 995
cor(Intercept,whether) 1344
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 31.68 1.61 28.41 34.87 1.00 2557 2982
whether 39.30 2.16 35.11 43.54 1.00 3175 2876
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 35.64 0.51 34.65 36.66 1.00 4587 2904
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Family: gaussian
Links: mu = identity; sigma = identity
Formula: rating ~ 1 + whether + how + (1 + whether + how | participant)
Data: df.data (Number of observations: 2624)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~participant (Number of levels: 41)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(Intercept) 5.27 2.01 1.24 9.36 1.00 1414
sd(whether) 13.10 2.41 8.76 18.05 1.00 1931
sd(how) 13.09 2.19 9.24 17.93 1.00 1141
cor(Intercept,whether) 0.29 0.32 -0.34 0.86 1.00 624
cor(Intercept,how) -0.18 0.32 -0.69 0.53 1.00 467
cor(whether,how) -0.63 0.16 -0.87 -0.25 1.00 757
Tail_ESS
sd(Intercept) 1535
sd(whether) 2994
sd(how) 1965
cor(Intercept,whether) 987
cor(Intercept,how) 640
cor(whether,how) 1864
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 10.29 1.56 7.30 13.41 1.00 3821 3177
whether 28.40 2.68 23.18 33.56 1.00 2444 3161
how 36.07 2.60 31.01 41.17 1.00 2066 2624
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 31.58 0.44 30.72 32.48 1.00 6371 3175
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Family: gaussian
Links: mu = identity; sigma = identity
Formula: rating ~ 1 + whether + how + sufficient + (1 + whether + how + sufficient | participant)
Data: df.data (Number of observations: 2624)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~participant (Number of levels: 41)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept) 6.27 2.08 1.93 10.42 1.00
sd(whether) 13.57 2.22 9.65 18.37 1.00
sd(how) 16.14 2.43 11.89 21.26 1.01
sd(sufficient) 14.49 2.84 9.35 20.52 1.00
cor(Intercept,whether) 0.12 0.27 -0.38 0.66 1.01
cor(Intercept,how) -0.16 0.26 -0.62 0.38 1.01
cor(whether,how) -0.71 0.13 -0.90 -0.42 1.01
cor(Intercept,sufficient) -0.01 0.28 -0.54 0.58 1.00
cor(whether,sufficient) 0.86 0.10 0.62 0.98 1.00
cor(how,sufficient) -0.68 0.15 -0.91 -0.34 1.00
Bulk_ESS Tail_ESS
sd(Intercept) 868 589
sd(whether) 1723 2357
sd(how) 1117 2189
sd(sufficient) 2004 2611
cor(Intercept,whether) 464 668
cor(Intercept,how) 438 794
cor(whether,how) 495 1618
cor(Intercept,sufficient) 565 845
cor(whether,sufficient) 1779 2769
cor(how,sufficient) 1543 2706
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 10.43 1.59 7.30 13.50 1.00 3044 3001
whether 27.93 2.60 22.80 33.09 1.00 1382 2074
how 26.47 2.91 20.83 32.42 1.00 1248 1933
sufficient 30.90 3.06 25.03 36.93 1.00 1643 2977
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 30.00 0.43 29.18 30.86 1.00 4986 3192
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
fit.brm_exp3_causal_w = add_criterion(fit.brm_exp3_causal_w,
criterion = c("loo", "waic"),
reloo = T)No problematic observations found. Returning the original 'loo' object.
fit.brm_exp3_causal_wh = add_criterion(fit.brm_exp3_causal_wh,
criterion = c("loo", "waic"),
reloo = T)No problematic observations found. Returning the original 'loo' object.
fit.brm_exp3_causal_whs = add_criterion(fit.brm_exp3_causal_whs,
criterion = c("loo", "waic"),
reloo = T)No problematic observations found. Returning the original 'loo' object.
loo_compare(fit.brm_exp3_causal_w,
fit.brm_exp3_causal_wh,
fit.brm_exp3_causal_whs)
# loo_compare(fit.brm_exp3_causal_w,
# fit.brm_exp3_causal_wh) elpd_diff se_diff
fit.brm_exp3_causal_whs 0.0 0.0
fit.brm_exp3_causal_wh -125.2 15.6
fit.brm_exp3_causal_w -426.2 31.1
# Regression based on features
df.regression.features = df.exp3.causal.long %>%
left_join(df.exp3.clipinfo, by = "clip") %>%
left_join(df.features %>%
mutate_at(.vars = vars(moving:exclusive), .funs = ~ scale(.)[,]),
by = c("clip", "ball")) %>%
clean_names() %>%
mutate(e_moving = 1 - e_moving)
# restrict the weights such that all predictors are positive
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))
fit.brm_exp3_causal_heuristic = brm(formula = rating ~ moving +
speed +
contact_e +
e_speed_diff +
e_direction_diff +
total_speed_diff +
total_direction_diff +
transfer +
e_moving +
exclusive + (1 | participant),
prior = priors,
data = df.regression.features %>%
select(-c(clip, ball, order, clipindex, contains("outcome"))),
save_all_pars = T,
cores = 2,
seed = 1,
file = "cache/brm_exp3_causal_heuristic")
fit.brm_exp3_causal_heuristic Family: gaussian
Links: mu = identity; sigma = identity
Formula: rating ~ moving + speed + contact_e + e_speed_diff + e_direction_diff + total_speed_diff + total_direction_diff + transfer + e_moving + exclusive + (1 | participant)
Data: df.regression.features %>% select(-c(clip, ball, o (Number of observations: 2624)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~participant (Number of levels: 41)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 8.53 1.22 6.44 11.18 1.00 1163 1936
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 49.79 1.46 46.91 52.74 1.00 1129
moving 0.22 0.22 0.01 0.81 1.00 3399
speed 2.06 0.85 0.44 3.67 1.00 2319
contact_e 0.38 0.35 0.01 1.29 1.00 3201
e_speed_diff 0.12 0.12 0.00 0.46 1.00 4435
e_direction_diff 1.04 0.74 0.04 2.71 1.00 2202
total_speed_diff 2.21 0.95 0.41 4.17 1.00 2208
total_direction_diff 3.99 0.91 2.14 5.68 1.00 2718
transfer 15.59 0.79 14.06 17.14 1.00 3555
e_moving 0.18 0.17 0.01 0.62 1.00 3296
exclusive 4.39 0.73 2.97 5.81 1.00 3526
Tail_ESS
Intercept 1854
moving 1547
speed 1266
contact_e 2020
e_speed_diff 2268
e_direction_diff 1685
total_speed_diff 1393
total_direction_diff 1742
transfer 2992
e_moving 1738
exclusive 2538
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 33.21 0.46 32.33 34.10 1.00 5313 2845
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
df.exp3.causal.regression = df.exp3.causal.means %>%
mutate(w = fit.brm_exp3_causal_w %>%
fitted(newdata = df.exp3.model.causal,
re_formula = NA) %>%
as_tibble() %>%
pull(Estimate),
wh = fit.brm_exp3_causal_wh %>%
fitted(newdata = df.exp3.model.causal,
re_formula = NA) %>%
as_tibble() %>%
pull(Estimate),
whs = fit.brm_exp3_causal_whs %>%
fitted(newdata = df.exp3.model.causal,
re_formula = NA) %>%
as_tibble() %>%
pull(Estimate),
heuristic = fit.brm_exp3_causal_heuristic %>%
fitted(newdata = df.exp3.model.causal %>%
left_join(df.features, by = c("clip", "ball")) %>%
clean_names(),
re_formula = NA) %>%
as_tibble() %>%
pull(Estimate))df.fit = df.exp3.causal.long %>%
left_join(df.exp3.model.causal,
by = c("clip", "ball"))
# priors
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))
# initial model fits (used for compilation)
fit.brm_exp3_causal_individual_baseline = brm(formula = rating ~ 1,
data = df.fit %>%
filter(participant == 1),
save_all_pars = T,
cores = 2,
seed = 1,
file = str_c("cache/brm_exp3_causal_individual_baseline"))
fit.brm_exp3_causal_individual_w = brm(formula = rating ~ 1 + whether,
data = df.fit %>%
filter(participant == 1),
prior = priors,
save_all_pars = T,
cores = 2,
seed = 1,
file = str_c("cache/brm_exp3_causal_individual_w"))
fit.brm_exp3_causal_individual_wh = brm(formula = rating ~ 1 + whether + how,
data = df.fit %>%
filter(participant == 1),
prior = priors,
save_all_pars = T,
cores = 2,
seed = 1,
file = str_c("cache/brm_exp3_causal_individual_wh"))
fit.brm_exp3_causal_individual_whs = brm(formula = rating ~ 1 + whether + how + sufficient,
data = df.fit %>%
filter(participant == 1),
prior = priors,
save_all_pars = T,
cores = 2,
seed = 1,
file = str_c("cache/brm_exp3_causal_individual_whs"))
# update model fits for different participants
df.exp3.causal.individual_fit = df.fit %>%
group_by(participant) %>%
nest() %>%
ungroup() %>%
# fit model for each participant
mutate(fit_baseline = map(.x = data,
.f = ~ update(fit.brm_exp3_causal_individual_baseline,
newdata = .x)),
fit_w = map(.x = data,
.f = ~ update(fit.brm_exp3_causal_individual_w,
newdata = .x)),
fit_wh = map(.x = data,
.f = ~ update(fit.brm_exp3_causal_individual_wh,
newdata = .x)),
fit_whs = map(.x = data,
.f = ~ update(fit.brm_exp3_causal_individual_whs,
newdata = .x))) %>%
mutate(fit_baseline = map(.x = fit_baseline,
~ add_criterion(.x, criterion = c("loo", "waic"))),
fit_w = map(.x = fit_w,
~ add_criterion(.x, criterion = c("loo", "waic"))),
fit_wh = map(.x = fit_wh,
~ add_criterion(.x, criterion = c("loo", "waic"))),
fit_whs = map(.x = fit_whs,
~ add_criterion(.x, criterion = c("loo", "waic"))),
r_w = map2_dbl(.x = data,
.y = fit_w,
.f = ~ cor(.x$rating, .y %>%
fitted() %>%
.[, 1])),
r_wh = map2_dbl(.x = data,
.y = fit_wh,
.f = ~ cor(.x$rating, .y %>%
fitted() %>%
.[, 1])),
r_whs = map2_dbl(.x = data,
.y = fit_whs,
.f = ~ cor(.x$rating, .y %>%
fitted() %>%
.[, 1])),
rmse_w = map2_dbl(.x = data,
.y = fit_w,
.f = ~ rmse(.x$rating, .y %>%
fitted() %>%
.[, 1])),
rmse_wh = map2_dbl(.x = data,
.y = fit_wh,
.f = ~ rmse(.x$rating, .y %>%
fitted() %>%
.[, 1])),
rmse_whs = map2_dbl(.x = data,
.y = fit_whs,
.f = ~ rmse(.x$rating, .y %>%
fitted() %>%
.[, 1])),
model_comparison = pmap(.l = list(baseline = fit_baseline,
w = fit_w,
wh = fit_wh,
whs = fit_whs),
.f = ~ loo_compare(..1, ..2, ..3, ..4)),
best_model = map_chr(.x = model_comparison,
.f = ~ rownames(.) %>% .[1]),
best_model = factor(best_model,
levels = c("..1", "..2", "..3", "..4"),
labels = c("baseline", "w", "wh", "whs")))
save(list = c("df.exp3.causal.individual_fit"),
file = "cache/exp3_causal_individual_fit.RData")load("cache/exp3_causal_individual_fit.RData")
# count how many participants are best fit by the different models
df.exp3.causal.individual_fit %>%
count(best_model)# A tibble: 2 x 2
best_model n
<fct> <int>
1 wh 2
2 whs 39
model_index = "whs"
# model predictions
model_prediction = list(fit.brm_exp3_causal_w,
fit.brm_exp3_causal_wh,
fit.brm_exp3_causal_whs) %>%
map(~ fitted(., df.exp3.causal.means %>%
left_join(df.exp3.model.causal,
by = c("clip", "ball")),
re_formula = NA)) %>%
reduce(rbind) %>%
as_tibble() %>%
mutate(ball = rep(c("A", "B"), n()/2),
clip = rep(rep(1:32, each = 2), 3),
model = rep(c("w", "wh", "whs"), each = 64))
df.plot = df.exp3.causal.long %>%
left_join(df.exp3.clipinfo %>%
select(clip, outcome_both),
by = c("clip")) %>%
left_join(model_prediction %>%
filter(model == model_index) %>%
select(model = Estimate,
clip,
ball),
by = c("clip", "ball")) %>%
gather("model","value",c(rating, model)) %>%
mutate(
colorindex = 1,
colorindex = ifelse(model == "rating" & outcome_both == 0, 1, colorindex),
colorindex = ifelse(model == "rating" & outcome_both == 1, 2, colorindex),
colorindex = ifelse(model == "model", 3, colorindex),
colorindex = as.factor(colorindex),
model = factor(model, levels = c("rating", "model"))
) %>%
arrange(participant, clip, ball)
df.text = df.plot %>%
expand(ball, clip, model) %>%
mutate(label = ifelse(model == "rating" & ball == "A", clip, NA),
y = 110,
x = 1.7,
colorindex = NA
)
actual_counterfactual_threeballs_plot(df.plot,"causal responsibility")
# ggsave("../../figures/plots/exp3_causal_bars.pdf",
# width = 12,
# height = 6)# \u2713 = check mark
# \u2717 = cross mark
check = "\u2713"
cross = "\u2717"
x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", cross),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))
df.plot = df.exp3.causal.long %>%
filter(clip %in% c(3, 7, 15, 16, 23)) %>%
group_by(clip, ball) %>%
summarise(mean = mean(rating),
low = smean.cl.boot(rating)[2],
high = smean.cl.boot(rating)[3]) %>%
left_join(df.exp3.causal.regression %>% select(clip, ball, w, wh, whs),
by = c("clip", "ball")) %>%
ungroup() %>%
gather("index", "value", c(mean, w, wh, whs)) %>%
mutate_at(.vars = vars(low, high), .funs = ~ ifelse(index == "mean", ., NA)) %>%
mutate(index = factor(index, levels = c("mean", "w", "wh", "whs")),
clip = factor(clip, levels = c(7, 23, 3, 15, 16),
labels = c("causal chain", "double prevention", "joint causation",
"overdetermination", "preemption")))
df.labels = df.plot %>%
distinct(clip, ball) %>%
arrange(clip, ball) %>%
mutate(labels = x_labels,
value = -10,
index = NA)
func_load_image = function(clip){
readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}
df.clips = df.plot %>%
distinct(clip) %>%
arrange(clip) %>%
mutate(number = c(7, 23, 3, 15, 16),
grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
index = NA,
value = NA,
ball = NA,
label = str_c("clip ", number))
ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)
ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)
p = ggplot(data = df.plot,
mapping = aes(x = ball,
y = value,
group = index,
fill = index)) +
geom_bar(stat = "identity", color = "black",
position = position_dodge(0.9),
width = 0.9) +
geom_errorbar(mapping = aes(ymin = low, ymax = high),
width = 0,
size = 1,
position = position_dodge(0.9)) +
annotation_custom(grob = ball_a, xmin = 0.5, xmax = 1.5, ymin = -25, ymax = -7) +
annotation_custom(grob = ball_b, xmin = 1.5, xmax = 2.5, ymin = -25, ymax = -7) +
geom_text(data = df.labels,
mapping = aes(label = labels,
y = -Inf),
vjust = 1.2,
family = "Arial Unicode MS",
size = 8) +
geom_custom(data = df.clips,
mapping = aes(data = grob, x = 1.5, y = Inf),
grob_fun = function(x) rasterGrob(x,
interpolate = T,
vjust = -0.25)) +
geom_text(data = df.clips,
mapping = aes(label = label,
y = Inf,
x = -Inf),
size = 9,
hjust = -0.2,
vjust = -4) +
facet_grid(~clip) +
scale_fill_grey(start = 1,
end = 0,
labels = c("mean rating ", expression(CSM[W] ~ " "),
expression(CSM[WH] ~ " "),
expression(CSM[WHS] ~ " "))) +
labs(y = "causal responsibility", fill = "") +
coord_cartesian(clip = "off") +
theme_bw() +
theme(legend.text.align = 0,
text = element_text(size = 20),
panel.grid = element_blank(),
legend.position = "bottom",
axis.title.x = element_blank(),
axis.text.x = element_blank(),
legend.text = element_text(size = 30),
strip.text = element_text(size = 30),
legend.spacing = unit(0.5, "cm"),
legend.background = element_rect(fill = "transparent"),
legend.margin = margin(t = 5, unit = "cm"),
plot.margin = margin(t = 8, l = 0.2, r = 0.2, b = 0.1, unit = "cm"))
# ggsave("../../figures/plots/exp3_selection_bars.pdf",
# width = 20,
# height = 10,
# device = cairo_pdf)model = "w"
model = "wh"
model = "whs"
model = "heuristic"
if(model == "heuristic"){
xlabel = "Heuristic"
}else{
xlabel = bquote(CSM[.(toupper(model))])
}
l.models = list(w = fit.brm_exp3_causal_w,
wh = fit.brm_exp3_causal_wh,
whs = fit.brm_exp3_causal_whs,
heuristic = fit.brm_exp3_causal_heuristic)
df.data = df.exp3.causal.means %>%
left_join(df.exp3.model.causal, by = c("clip", "ball")) %>%
left_join(df.regression.features, by = c("clip", "ball"))
df.plot = l.models[[model]] %>%
fitted(newdata = df.data,
re_formula = NA) %>%
as_tibble() %>%
clean_names() %>%
bind_cols(df.data)
ggplot(data = df.plot,
mapping = aes(x = estimate,
y = rating_mean)) +
geom_abline(intercept = 0,
slope = 1,
linetype = 2) +
geom_smooth(aes(y = estimate,
ymin = q2_5,
ymax = q97_5),
stat = "identity",
color = "black") +
geom_linerange(size = 0.75,
mapping = aes(ymin = rating_low,
ymax = rating_high),
color = "gray50") +
geom_point(size = 2) +
scale_color_manual(values = c("black", "#e41a1c"),
guide = F) +
coord_fixed(xlim = c(0, 100),
ylim = c(0, 100)) +
labs(x = xlabel,
y = "responsibility judgment") +
annotate(geom = "text",
label = str_c(
"r = ", cor(df.plot$estimate, df.plot$rating_mean) %>% round(2) %>% as.character() %>% str_sub(start = 2),
"\nRMSE = ", rmse(df.plot$estimate, df.plot$rating_mean) %>% round(2)),
x = 5,
y = 95,
size = 6,
hjust = 0) +
scale_x_continuous(breaks = seq(0, 100, 25)) +
scale_y_continuous(breaks = seq(0, 100, 25)) +
theme(text = element_text(size = 20))
# ggsave(str_c("../../figures/plots/exp3_", model, "_scatter.pdf"),
# width = 5,
# height = 5)df.plot = df.exp3.causal.individual_fit %>%
select(participant, fit_whs) %>%
mutate(estimates = map(fit_whs, tidy)) %>%
select(participant, estimates) %>%
unnest(estimates) %>%
filter(str_detect(term, "b_"),
!str_detect(term, "Intercept")) %>%
mutate(term = str_remove(term, "b_")) %>%
select(participant, term, estimate) %>%
spread(term, estimate) %>%
mutate_at(vars(how:whether), .funs = list(norm = ~ . / (how + whether + sufficient))) %>%
mutate(color = 0,
color = ifelse(how_norm == max(how_norm), 1, color),
color = ifelse(whether_norm == max(whether_norm), 2, color),
color = ifelse(sufficient_norm == max(sufficient_norm), 3, color),
color = factor(color))
ggplot(data = df.plot,
mapping = aes(x = how,
y = sufficient,
z = whether,
color = color)) +
geom_point(alpha = 0.7,
size = 2,
show.legend = F) +
scale_color_manual(values = c("black", "red", "green", "blue")) +
coord_tern() +
theme_bw() +
theme_showarrows() +
theme(text = element_text(size = 20),
tern.axis.title.T = element_blank(),
tern.axis.title.L = element_blank(),
tern.axis.title.R = element_blank())
# ggsave(str_c("../../figures/plots/exp3_individual_regression_ternary_plot_scaled.pdf"),
# width = 5,
# height = 5)# \u2713 = check mark
# \u2717 = cross mark
check = "\u2713"
cross = "\u2717"
x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", cross),
str_c("\nW ", check, "\nH ", check, "\nS ", cross),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", check, "\nS ", check),
str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))
df.plot = df.exp3.causal.long %>%
filter(clip %in% c(7, 23, 3, 15, 16)) %>%
mutate(clip = factor(clip, levels = c(7, 23, 3, 15, 16),
labels = c("causal chain", "double prevention", "joint causation",
"overdetermination", "preemption")))
df.cluster = df.exp3.causal.individual_fit %>%
select(participant, fit_whs) %>%
mutate(estimates = map(fit_whs, tidy)) %>%
select(participant, estimates) %>%
unnest(estimates) %>%
filter(str_detect(term, "b_"),
!str_detect(term, "Intercept")) %>%
mutate(term = str_remove(term, "b_")) %>%
select(participant, term, estimate) %>%
spread(term, estimate) %>%
mutate(group = ifelse(how > whether, "how", "whether"))
df.plot = df.plot %>%
left_join(df.cluster %>%
select(participant, group),
by = "participant")
df.labels = df.plot %>%
distinct(clip, ball) %>%
arrange(clip, ball) %>%
mutate(labels = x_labels,
rating = -10,
group = NA,
participant = NA)
func_load_image = function(clip){
readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}
df.clips = df.plot %>%
distinct(clip) %>%
arrange(clip) %>%
mutate(number = c(7, 23, 3, 15, 16),
grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
group = NA,
participant = NA,
ball = NA,
label = str_c("clip ", number))
ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)
ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)
# ggplot(df.plot,aes(x=ball,y=rating,group=id,color = best))+
p = ggplot(data = df.plot,
mapping = aes(x = ball,
y = rating,
group = participant,
shape = group)) +
geom_line(mapping = aes(linetype = "individual",
color = group),
size = 1,
alpha = 0.3) +
geom_point(mapping = aes(color = group),
size = 1,
alpha = 0.3) +
stat_summary(mapping = aes(group = group,
color = group,
linetype = "mean"),
fun.y = "mean",
geom = "line",
size = 1.5) +
stat_summary(data = df.plot %>% filter(group == "how"),
mapping = aes(group = group,
color = group),
fun.data = "mean_cl_boot",
geom = "pointrange",
size = 1,
shape = 19) +
stat_summary(data = df.plot %>% filter(group == "whether"),
mapping = aes(group = group,
color = group),
fun.data = "mean_cl_boot",
geom = "pointrange",
size = 1,
shape = 19) +
annotation_custom(grob = ball_a, xmin = 0.5, xmax = 1.5, ymin = -30, ymax = -10) +
annotation_custom(grob = ball_b, xmin = 1.5, xmax = 2.5, ymin = -30, ymax = -10) +
geom_text(data = df.labels,
mapping = aes(label = labels,
y = -Inf),
vjust = 1.2,
family = "Arial Unicode MS",
size = 8) +
geom_custom(data = df.clips,
mapping = aes(data = grob, x = 1.5, y = Inf),
grob_fun = function(x) rasterGrob(x,
interpolate = T,
vjust = -0.25)) +
geom_text(data = df.clips,
mapping = aes(label = label,
y = Inf,
x = -Inf),
size = 9,
hjust = -0.2,
vjust = -4) +
facet_wrap(~clip,
ncol = 8) +
coord_cartesian(xlim = c(0.9, 2.1),
ylim = c(-5, 105),
expand = F,
clip = "off") +
scale_y_continuous(breaks = seq(0, 100, 25),
labels = seq(0, 100, 25)) +
scale_color_brewer(palette = "Set1",
guide = "none") +
labs(y = "causal responsibility rating",
linetype = "legend",
shape = "") +
theme_bw() +
guides(linetype = guide_legend(override.aes = list(alpha = c(0.3, 1)),
keywidth = unit(1.2, "cm")),
shape = guide_legend(override.aes = list(color = c("#e41a1c", "#377eb8"),
shape = c(19, 19),
alpha = c(1, 1),
size = c(5, 5)))) +
theme(panel.grid = element_blank(),
text = element_text(size = 20),
legend.position = c(0.505, 0.23),
axis.title.x = element_blank(),
legend.text = element_text(size = 20),
legend.title = element_text(size = 24),
legend.background = element_rect(fill = "transparent"),
strip.text = element_text(size = 30),
axis.text.x = element_blank(),
legend.key = element_rect(fill = "transparent"),
legend.box = "horizontal",
legend.spacing.x = unit(0.1, "cm"),
panel.spacing.x = unit(1, "cm"),
plot.margin = margin(b = 5, l = 0.2, r = 0.2, t = 7.5, unit = "cm"))
# ggsave(str_c("../../figures/plots/exp3_individual_variance_selection_lines_clustered.pdf"),
# plot = p,
# width = 20,
# height = 8.5,
# device = cairo_pdf)df.exp3.model %>%
filter(noise == unique(df.exp3.model.counterfactual$noise)) %>%
gather("index", "value", A_difference:B_robust) %>%
separate(index, into = c("ball", "aspect")) %>%
spread(aspect, value) %>%
select(difference,
whether,
how,
sufficient,
robust) %>%
correlate() %>%
shave() %>%
fashion() %>%
xtable() %>%
print(include.rownames = F,
booktabs = T)
Correlation method: 'pearson'
Missing treated using: 'pairwise.complete.obs'
% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:06:10 2019
\begin{table}[ht]
\centering
\begin{tabular}{llllll}
\toprule
rowname & difference & whether & how & sufficient & robust \\
\midrule
difference & & & & & \\
whether & .48 & & & & \\
how & .76 & .24 & & & \\
sufficient & .18 & .05 & .34 & & \\
robust & .42 & .94 & .21 & .13 & \\
\bottomrule
\end{tabular}
\end{table}
fit.brm_exp3_causal_w %>%
tidy() %>%
filter(str_detect(term, "b_")) %>%
mutate(model = "CSM_w") %>%
bind_rows(fit.brm_exp3_causal_wh %>%
tidy() %>%
filter(str_detect(term, "b_")) %>%
mutate(model = "CSM_wh")) %>%
bind_rows(fit.brm_exp3_causal_whs %>%
tidy() %>%
filter(str_detect(term, "b_")) %>%
mutate(model = "CSM_whs")) %>%
mutate(term = tolower(term),
term = str_remove_all(term, "b_")) %>%
rename(`lower 95% HDI` = lower,
`upper 95% HDI` = upper) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
select(model, everything()) %>%
xtable() %>%
print(include.rownames = F,
booktabs = T)% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:06:11 2019
\begin{table}[ht]
\centering
\begin{tabular}{llrrrr}
\toprule
model & term & estimate & std.error & lower 95\% HDI & upper 95\% HDI \\
\midrule
CSM\_w & intercept & 31.68 & 1.61 & 29.02 & 34.37 \\
CSM\_w & whether & 39.30 & 2.16 & 35.76 & 42.82 \\
CSM\_wh & intercept & 10.29 & 1.56 & 7.80 & 12.83 \\
CSM\_wh & whether & 28.40 & 2.68 & 23.96 & 32.70 \\
CSM\_wh & how & 36.07 & 2.60 & 31.80 & 40.27 \\
CSM\_whs & intercept & 10.43 & 1.59 & 7.82 & 13.01 \\
CSM\_whs & whether & 27.93 & 2.60 & 23.67 & 32.18 \\
CSM\_whs & how & 26.47 & 2.91 & 21.76 & 31.35 \\
CSM\_whs & sufficient & 30.90 & 3.06 & 25.77 & 35.96 \\
\bottomrule
\end{tabular}
\end{table}
df.exp3.causal.individual_fit %>%
select_if(~ (is.numeric(.) | is.factor(.))) %>%
mutate_at(.vars = vars(contains("_w")), .funs = ~ round(., 2)) %>%
select(participant, everything(), best_model) %>%
xtable() %>%
print(include.rownames = F,
booktabs = T)% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:06:12 2019
\begin{table}[ht]
\centering
\begin{tabular}{rrrrrrrl}
\toprule
participant & r\_w & r\_wh & r\_whs & rmse\_w & rmse\_wh & rmse\_whs & best\_model \\
\midrule
1 & 0.29 & 0.38 & 0.48 & 30.00 & 29.03 & 27.74 & whs \\
2 & 0.20 & 0.77 & 0.77 & 39.32 & 27.94 & 27.55 & whs \\
3 & 0.37 & 0.78 & 0.79 & 38.90 & 28.00 & 26.95 & whs \\
4 & 0.41 & 0.55 & 0.63 & 43.62 & 40.60 & 38.31 & whs \\
5 & 0.45 & 0.53 & 0.54 & 39.23 & 37.32 & 36.81 & whs \\
6 & 0.18 & 0.54 & 0.54 & 40.84 & 36.32 & 36.07 & whs \\
7 & 0.49 & 0.61 & 0.64 & 32.96 & 29.93 & 29.15 & whs \\
8 & 0.39 & 0.63 & 0.68 & 38.37 & 33.26 & 31.64 & whs \\
9 & 0.40 & 0.72 & 0.76 & 35.86 & 28.34 & 26.39 & whs \\
10 & 0.28 & 0.52 & 0.53 & 29.13 & 26.28 & 25.85 & whs \\
11 & 0.51 & 0.56 & 0.60 & 33.66 & 32.19 & 31.17 & whs \\
12 & 0.56 & 0.57 & 0.70 & 32.78 & 32.08 & 28.84 & whs \\
13 & 0.56 & 0.71 & 0.74 & 29.83 & 25.40 & 24.21 & whs \\
14 & 0.43 & 0.47 & 0.50 & 33.72 & 32.82 & 32.23 & whs \\
15 & 0.19 & 0.35 & 0.37 & 39.77 & 38.29 & 37.94 & whs \\
16 & 0.62 & 0.67 & 0.76 & 40.49 & 37.98 & 34.39 & whs \\
17 & 0.42 & 0.66 & 0.71 & 28.57 & 23.82 & 22.42 & whs \\
18 & 0.31 & 0.35 & 0.39 & 37.78 & 37.09 & 36.57 & whs \\
19 & 0.37 & 0.60 & 0.62 & 33.41 & 29.23 & 28.62 & whs \\
20 & 0.57 & 0.61 & 0.68 & 36.47 & 35.00 & 32.71 & whs \\
21 & 0.43 & 0.50 & 0.57 & 30.95 & 29.61 & 28.32 & whs \\
22 & 0.51 & 0.67 & 0.71 & 33.72 & 29.36 & 28.00 & whs \\
23 & 0.33 & 0.45 & 0.52 & 38.76 & 37.00 & 35.57 & whs \\
24 & 0.46 & 0.63 & 0.69 & 34.01 & 29.96 & 28.22 & whs \\
25 & 0.61 & 0.66 & 0.70 & 30.60 & 28.91 & 27.55 & whs \\
26 & 0.57 & 0.71 & 0.76 & 40.43 & 35.03 & 32.80 & whs \\
27 & 0.46 & 0.52 & 0.66 & 43.12 & 41.49 & 38.12 & whs \\
28 & 0.30 & 0.45 & 0.46 & 29.33 & 27.62 & 27.38 & whs \\
29 & 0.48 & 0.60 & 0.65 & 38.93 & 36.03 & 34.48 & whs \\
30 & 0.18 & 0.70 & 0.70 & 33.06 & 24.95 & 24.71 & whs \\
31 & 0.48 & 0.59 & 0.69 & 37.54 & 34.72 & 31.76 & whs \\
32 & 0.37 & 0.61 & 0.70 & 40.62 & 35.94 & 32.92 & whs \\
33 & 0.50 & 0.56 & 0.64 & 38.64 & 36.79 & 34.85 & whs \\
34 & 0.24 & 0.51 & 0.52 & 44.50 & 40.86 & 40.19 & whs \\
35 & 0.15 & 0.32 & 0.30 & 34.56 & 33.36 & 33.37 & wh \\
36 & 0.19 & 0.84 & 0.82 & 45.38 & 29.58 & 29.70 & wh \\
37 & 0.55 & 0.63 & 0.67 & 36.04 & 33.36 & 32.11 & whs \\
38 & 0.52 & 0.61 & 0.67 & 36.03 & 33.39 & 31.53 & whs \\
39 & 0.47 & 0.77 & 0.77 & 29.61 & 21.75 & 21.58 & whs \\
40 & 0.46 & 0.78 & 0.79 & 32.41 & 23.60 & 22.69 & whs \\
41 & 0.32 & 0.50 & 0.58 & 32.90 & 30.38 & 28.69 & whs \\
\bottomrule
\end{tabular}
\end{table}
df.table = df.exp3.model.causal %>%
left_join(df.exp3.causal.regression) %>%
filter(clip %in% c(7, 23, 16, 3, 15)) %>%
mutate_at(.vars = vars(difference, robust, sufficient, whether, heuristic), .funs = ~ round(., 2)) %>%
select(clip, ball, difference, whether, how, sufficient, robust) %>%
mutate(clip = factor(clip, levels = c(7, 23, 16, 3, 15))) %>%
mutate_at(.vars = vars(difference, whether, how, sufficient, robust), .funs = ~ ifelse(. < 0.5, paste0("xmark (", ., ")"), paste0("cmark (", ., ")"))) %>%
arrange(clip)Joining, by = c("clip", "ball")
% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:06:12 2019
\begin{table}[ht]
\centering
\begin{tabular}{lllllll}
\hline
clip & ball & difference & whether & how & sufficient & robust \\
\hline
7 & A & cmark (1) & xmark (0.2) & cmark (1) & xmark (0) & xmark (0.18) \\
7 & B & cmark (1) & cmark (1) & cmark (1) & cmark (0.81) & cmark (0.6) \\
23 & A & xmark (0.05) & xmark (0) & xmark (0) & xmark (0) & xmark (0) \\
23 & B & cmark (0.96) & cmark (0.89) & xmark (0) & xmark (0) & cmark (0.8) \\
16 & A & cmark (1) & xmark (0.19) & cmark (1) & cmark (1) & xmark (0.29) \\
16 & B & xmark (0) & xmark (0) & xmark (0) & xmark (0) & xmark (0) \\
3 & A & cmark (1) & cmark (0.88) & cmark (1) & xmark (0.13) & cmark (0.74) \\
3 & B & cmark (1) & cmark (0.87) & cmark (1) & xmark (0.12) & cmark (0.74) \\
15 & A & cmark (1) & xmark (0) & cmark (1) & cmark (1) & xmark (0.08) \\
15 & B & cmark (1) & xmark (0) & cmark (1) & cmark (1) & xmark (0.06) \\
\hline
\end{tabular}
\end{table}
df.table = df.exp3.causal.regression %>%
left_join(df.exp3.model.causal,
by = c("clip", "ball")) %>%
left_join(df.exp3.clipinfo %>%
select(-clipindex),
by = c("clip")) %>%
mutate_at(.vars = vars(difference, whether, how, sufficient, robust),
.funs = ~ . * 100) %>%
mutate_at(.vars = vars(-ball), .funs = ~ round(., 0)) %>%
select(clip, ball, contains("outcome"), difference, whether, how, sufficient, robust, w, wh, whs, heuristic, rating = rating_mean)
xtable(df.table, digits = 0) %>%
print(include.rownames = FALSE)% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:06:12 2019
\begin{table}[ht]
\centering
\begin{tabular}{rlrrrrrrrrrrrrrr}
\hline
clip & ball & outcome\_both & outcome\_a & outcome\_b & outcome\_none & difference & whether & how & sufficient & robust & w & wh & whs & heuristic & rating \\
\hline
1 & A & 0 & 0 & 0 & 0 & 100 & 41 & 100 & 20 & 40 & 48 & 58 & 54 & 75 & 42 \\
1 & B & 0 & 0 & 0 & 0 & 100 & 19 & 100 & 16 & 7 & 39 & 52 & 47 & 72 & 37 \\
2 & A & 0 & 0 & 0 & 0 & 58 & 11 & 0 & 0 & 12 & 36 & 13 & 14 & 64 & 21 \\
2 & B & 0 & 0 & 0 & 0 & 14 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 64 & 19 \\
3 & A & 1 & 0 & 0 & 0 & 100 & 88 & 100 & 13 & 74 & 66 & 71 & 66 & 92 & 76 \\
3 & B & 1 & 0 & 0 & 0 & 100 & 87 & 100 & 12 & 74 & 66 & 71 & 65 & 92 & 75 \\
4 & A & 1 & 0 & 0 & 0 & 100 & 76 & 100 & 0 & 70 & 62 & 68 & 58 & 76 & 63 \\
4 & B & 1 & 0 & 0 & 0 & 100 & 99 & 100 & 16 & 60 & 71 & 74 & 69 & 72 & 78 \\
5 & A & 0 & 0 & 1 & 0 & 100 & 97 & 100 & 0 & 49 & 70 & 74 & 64 & 68 & 71 \\
5 & B & 0 & 0 & 1 & 0 & 100 & 0 & 100 & 0 & 0 & 32 & 46 & 37 & 88 & 22 \\
6 & A & 0 & 0 & 1 & 0 & 100 & 67 & 100 & 12 & 38 & 58 & 66 & 59 & 71 & 73 \\
6 & B & 0 & 0 & 1 & 0 & 100 & 12 & 100 & 4 & 12 & 36 & 50 & 42 & 72 & 22 \\
7 & A & 1 & 0 & 1 & 0 & 100 & 20 & 100 & 0 & 18 & 40 & 52 & 43 & 86 & 59 \\
7 & B & 1 & 0 & 1 & 0 & 100 & 100 & 100 & 81 & 60 & 71 & 75 & 90 & 85 & 79 \\
8 & A & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 65 & 7 \\
8 & B & 1 & 0 & 1 & 0 & 100 & 100 & 100 & 100 & 100 & 71 & 75 & 96 & 98 & 92 \\
9 & A & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 53 & 8 \\
9 & B & 0 & 1 & 0 & 0 & 100 & 100 & 100 & 0 & 100 & 71 & 75 & 65 & 90 & 90 \\
10 & A & 0 & 1 & 0 & 0 & 82 & 13 & 0 & 0 & 17 & 37 & 14 & 14 & 55 & 23 \\
10 & B & 0 & 1 & 0 & 0 & 100 & 88 & 0 & 0 & 68 & 66 & 35 & 35 & 59 & 55 \\
11 & A & 1 & 1 & 0 & 0 & 100 & 74 & 100 & 76 & 70 & 61 & 67 & 81 & 82 & 93 \\
11 & B & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 54 & 4 \\
12 & A & 1 & 1 & 0 & 0 & 100 & 84 & 100 & 75 & 81 & 65 & 70 & 84 & 75 & 77 \\
12 & B & 1 & 1 & 0 & 0 & 100 & 0 & 100 & 9 & 25 & 32 & 46 & 40 & 72 & 37 \\
13 & A & 0 & 1 & 1 & 0 & 71 & 44 & 0 & 0 & 40 & 49 & 23 & 23 & 54 & 8 \\
13 & B & 0 & 1 & 1 & 0 & 81 & 51 & 0 & 0 & 47 & 52 & 25 & 25 & 59 & 64 \\
14 & A & 0 & 1 & 1 & 0 & 100 & 97 & 0 & 0 & 60 & 70 & 38 & 38 & 59 & 22 \\
14 & B & 0 & 1 & 1 & 0 & 96 & 86 & 0 & 0 & 58 & 66 & 35 & 34 & 58 & 18 \\
15 & A & 1 & 1 & 1 & 0 & 100 & 0 & 100 & 100 & 8 & 32 & 46 & 68 & 91 & 76 \\
15 & B & 1 & 1 & 1 & 0 & 100 & 0 & 100 & 100 & 6 & 32 & 46 & 68 & 91 & 76 \\
16 & A & 1 & 1 & 1 & 0 & 100 & 19 & 100 & 100 & 29 & 39 & 52 & 73 & 94 & 92 \\
16 & B & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 53 & 4 \\
17 & A & 0 & 0 & 0 & 1 & 100 & 11 & 100 & 53 & 14 & 36 & 49 & 56 & 86 & 69 \\
17 & B & 0 & 0 & 0 & 1 & 100 & 0 & 100 & 45 & 18 & 32 & 46 & 51 & 84 & 46 \\
18 & A & 0 & 0 & 0 & 1 & 100 & 6 & 100 & 52 & 12 & 34 & 48 & 54 & 74 & 63 \\
18 & B & 0 & 0 & 0 & 1 & 100 & 4 & 100 & 50 & 6 & 33 & 48 & 54 & 74 & 66 \\
19 & A & 1 & 0 & 0 & 1 & 100 & 78 & 100 & 3 & 71 & 62 & 69 & 60 & 73 & 53 \\
19 & B & 1 & 0 & 0 & 1 & 100 & 80 & 100 & 4 & 71 & 63 & 69 & 60 & 73 & 49 \\
20 & A & 1 & 0 & 0 & 1 & 100 & 92 & 100 & 3 & 72 & 68 & 73 & 64 & 75 & 41 \\
20 & B & 1 & 0 & 0 & 1 & 100 & 94 & 100 & 3 & 59 & 69 & 73 & 64 & 74 & 71 \\
21 & A & 0 & 0 & 1 & 1 & 100 & 56 & 100 & 53 & 53 & 54 & 62 & 69 & 76 & 80 \\
21 & B & 0 & 0 & 1 & 1 & 100 & 4 & 100 & 24 & 4 & 33 & 47 & 45 & 78 & 18 \\
22 & A & 0 & 0 & 1 & 1 & 100 & 100 & 100 & 96 & 84 & 71 & 75 & 95 & 69 & 60 \\
22 & B & 0 & 0 & 1 & 1 & 100 & 4 & 100 & 0 & 14 & 33 & 48 & 38 & 71 & 42 \\
23 & A & 1 & 0 & 1 & 1 & 5 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 54 & 3 \\
23 & B & 1 & 0 & 1 & 1 & 96 & 89 & 0 & 0 & 80 & 67 & 35 & 35 & 60 & 39 \\
24 & A & 1 & 0 & 1 & 1 & 100 & 61 & 100 & 1 & 54 & 56 & 64 & 54 & 75 & 44 \\
24 & B & 1 & 0 & 1 & 1 & 100 & 98 & 100 & 21 & 84 & 70 & 74 & 71 & 72 & 73 \\
25 & A & 0 & 1 & 0 & 1 & 100 & 18 & 100 & 17 & 19 & 39 & 52 & 47 & 80 & 43 \\
25 & B & 0 & 1 & 0 & 1 & 100 & 78 & 100 & 64 & 69 & 62 & 68 & 78 & 75 & 73 \\
26 & A & 0 & 1 & 0 & 1 & 100 & 10 & 100 & 3 & 12 & 36 & 49 & 41 & 78 & 39 \\
26 & B & 0 & 1 & 0 & 1 & 100 & 96 & 100 & 42 & 57 & 69 & 74 & 77 & 68 & 69 \\
27 & A & 1 & 1 & 0 & 1 & 100 & 100 & 100 & 47 & 98 & 71 & 75 & 79 & 78 & 80 \\
27 & B & 1 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 55 & 6 \\
28 & A & 1 & 1 & 0 & 1 & 100 & 94 & 100 & 14 & 84 & 69 & 73 & 68 & 92 & 89 \\
28 & B & 1 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 53 & 5 \\
29 & A & 0 & 1 & 1 & 1 & 100 & 64 & 100 & 21 & 49 & 57 & 65 & 61 & 78 & 47 \\
29 & B & 0 & 1 & 1 & 1 & 100 & 70 & 100 & 23 & 43 & 59 & 66 & 64 & 73 & 67 \\
30 & A & 0 & 1 & 1 & 1 & 100 & 65 & 100 & 28 & 55 & 57 & 65 & 64 & 83 & 58 \\
30 & B & 0 & 1 & 1 & 1 & 100 & 55 & 100 & 24 & 45 & 53 & 62 & 60 & 78 & 56 \\
31 & A & 1 & 1 & 1 & 1 & 100 & 0 & 100 & 1 & 0 & 32 & 46 & 37 & 75 & 44 \\
31 & B & 1 & 1 & 1 & 1 & 100 & 0 & 100 & 1 & 0 & 32 & 47 & 37 & 71 & 46 \\
32 & A & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 54 & 5 \\
32 & B & 1 & 1 & 1 & 1 & 100 & 66 & 100 & 52 & 65 & 58 & 65 & 71 & 77 & 71 \\
\hline
\end{tabular}
\end{table}
fit.brm_exp3_causal_heuristic %>%
tidy() %>%
filter(str_detect(term, "b_")) %>%
mutate(term = str_remove(term, "b_"),
term = tolower(term)) %>%
mutate_if(is.numeric, ~ round(., 2)) %>%
rename(`lower 95% HDI` = lower,
`upper 95% HDI` = upper) %>%
xtable() %>%
print(include.rownames = F,
booktabs = T)% latex table generated in R 3.6.1 by xtable 1.8-4 package
% Tue Nov 12 15:06:12 2019
\begin{table}[ht]
\centering
\begin{tabular}{lrrrr}
\toprule
term & estimate & std.error & lower 95\% HDI & upper 95\% HDI \\
\midrule
intercept & 49.79 & 1.46 & 47.39 & 52.23 \\
moving & 0.22 & 0.22 & 0.01 & 0.66 \\
speed & 2.06 & 0.85 & 0.64 & 3.43 \\
contact\_e & 0.38 & 0.35 & 0.02 & 1.07 \\
e\_speed\_diff & 0.12 & 0.12 & 0.01 & 0.37 \\
e\_direction\_diff & 1.04 & 0.74 & 0.08 & 2.43 \\
total\_speed\_diff & 2.21 & 0.95 & 0.67 & 3.80 \\
total\_direction\_diff & 3.99 & 0.91 & 2.45 & 5.45 \\
transfer & 15.59 & 0.79 & 14.30 & 16.90 \\
e\_moving & 0.18 & 0.17 & 0.01 & 0.51 \\
exclusive & 4.39 & 0.73 & 3.16 & 5.56 \\
\bottomrule
\end{tabular}
\end{table}
R version 3.6.1 (2019-07-05)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] grid stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] forcats_0.4.0 stringr_1.4.0 dplyr_0.8.3
[4] purrr_0.3.2 readr_1.3.1 tidyr_1.0.0
[7] tibble_2.1.3 tidyverse_1.2.1 egg_0.4.5
[10] gridExtra_2.3 png_0.1-7 tidybayes_1.1.0
[13] Hmisc_4.2-0 Formula_1.2-3 survival_2.44-1.1
[16] lattice_0.20-38 ggtern_3.1.0 ggplot2_3.2.1
[19] janitor_1.2.0 jsonlite_1.6 xtable_1.8-4
[22] corrr_0.4.0 brms_2.10.0 Rcpp_1.0.2
[25] broom_0.5.2 lme4_1.1-21 Matrix_1.2-17
[28] knitr_1.25
loaded via a namespace (and not attached):
[1] readxl_1.3.1 backports_1.1.4
[3] plyr_1.8.4 igraph_1.2.4.1
[5] lazyeval_0.2.2 splines_3.6.1
[7] svUnit_0.7-12 crosstalk_1.0.0
[9] rstantools_2.0.0 inline_0.3.15
[11] digest_0.6.21 htmltools_0.3.6
[13] rsconnect_0.8.15 fansi_0.4.0
[15] magrittr_1.5 checkmate_1.9.4
[17] cluster_2.1.0 modelr_0.1.5
[19] bayesm_3.1-3 matrixStats_0.55.0
[21] xts_0.11-2 prettyunits_1.0.2
[23] colorspace_1.4-1 rvest_0.3.4
[25] haven_2.1.1 xfun_0.9
[27] callr_3.3.2 crayon_1.3.4
[29] zeallot_0.1.0 zoo_1.8-6
[31] glue_1.3.1 gtable_0.3.0
[33] compositions_1.40-2 pkgbuild_1.0.5
[35] rstan_2.19.2 DEoptimR_1.0-8
[37] abind_1.4-5 scales_1.0.0
[39] miniUI_0.1.1.1 htmlTable_1.13.2
[41] HDInterval_0.2.0 ggstance_0.3.3
[43] foreign_0.8-72 latex2exp_0.4.0
[45] stats4_3.6.1 StanHeaders_2.19.0
[47] DT_0.9 htmlwidgets_1.3
[49] httr_1.4.1 threejs_0.3.1
[51] arrayhelpers_1.0-20160527 RColorBrewer_1.1-2
[53] ellipsis_0.3.0 acepack_1.4.1
[55] pkgconfig_2.0.3 loo_2.1.0
[57] nnet_7.3-12 utf8_1.1.4
[59] labeling_0.3 tidyselect_0.2.5
[61] rlang_0.4.0 reshape2_1.4.3
[63] later_0.8.0 cellranger_1.1.0
[65] munsell_0.5.0 tools_3.6.1
[67] cli_1.1.0 generics_0.0.2
[69] ggridges_0.5.1 evaluate_0.14
[71] yaml_2.2.0 processx_3.4.1
[73] robustbase_0.93-5 nlme_3.1-141
[75] mime_0.7 xml2_1.2.2
[77] compiler_3.6.1 bayesplot_1.7.0
[79] shinythemes_1.1.2 rstudioapi_0.10
[81] stringi_1.4.3 ps_1.3.0
[83] Brobdingnag_1.2-6 nloptr_1.2.1
[85] markdown_1.1 shinyjs_1.0
[87] tensorA_0.36.1 vctrs_0.2.0
[89] pillar_1.4.2 lifecycle_0.1.0
[91] bridgesampling_0.7-2 data.table_1.12.2
[93] httpuv_1.5.2 R6_2.4.0
[95] latticeExtra_0.6-28 bookdown_0.13
[97] promises_1.0.1 boot_1.3-23
[99] energy_1.7-6 colourpicker_1.0
[101] MASS_7.3-51.4 gtools_3.8.1
[103] assertthat_0.2.1 proto_1.0.0
[105] withr_2.1.2 shinystan_2.5.0
[107] parallel_3.6.1 hms_0.5.1
[109] rpart_4.1-15 coda_0.19-3
[111] minqa_1.2.4 snakecase_0.11.0
[113] rmarkdown_1.15 shiny_1.3.2
[115] lubridate_1.7.4 base64enc_0.1-3
[117] dygraphs_1.1.1.6